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Abstract 

This  study  uses  a  linear  model  of  an  Integrated  Power  and  Attitude  Control 
System  (IPACS)  to  investigate  the  vibration  interaction  between  multiple  flywheels. 
An  easily  extendable  Matlab®  script  is  created  for  the  analysis  of  flywheel  vibra¬ 
tions.  This  script  is  used  to  build  a  vibration  model  consisting  of  two  active  magnetic 
bearing  flywheels  mounted  on  a  support  structure.  The  flywheels  are  rotated  at  vary¬ 
ing  speeds,  with  an  imbalance-induced  centripetal  force  in  one  or  both  wheels  causing 
vibrations  in  both  wheels.  Flywheel  and  system  responses  are  examined  for  low  fre¬ 
quency  vibrations  which  would  cause  undesirable  excitation  to  a  satellite  using  IPACS, 
with  a  specihc  focus  on  the  beat  phenomenon  and  extra-synchronous  vibration.  Extra- 
synchronous  resonant  vibration  between  multiple  rotors  is  shown  to  exist  in  an  ideal 
undamped  conhguration  but  even  a  very  small  realistic  amount  of  damping  is  enough 
to  mitigate  the  effect  enough  that  it  is  of  less  concern  than  individual  rotor  vibration 
inputs.  Extra-synchronous  resonant  vibration  is  thus  shown  to  have  a  minimal  effect 
on  satellite  IPACS  operation. 
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Vibration  Interaction  in  a 
Multiple  Flywheel  System 

I.  Introduction 

Advanced  flywheels  are  an  exciting  technology  with  the  potential  to  greatly  im¬ 
prove  performance  for  satellite  energy  storage  systems.  They  have  been  investigated 
for  use  in  space  since  the  early  1960 ’s,  and  supporting  technologies  are  Anally  begin¬ 
ning  to  mature  to  the  point  that  they  may  soon  be  feasible.  Unfortunately,  despite 
much  historic  optimism,  there  are  still  unsolved  and  unstudied  problems  with  their 
operation  and  implementation.  This  thesis  investigates  two  areas  of  potential  concern 
to  see  whether  they  pose  any  particular  challenges  to  advanced  flywheel  operation  in 
space:  the  beat  phenomenon  and  extra-synchronous  whirl  excitation  caused  by  inter¬ 
actions  between  multiple  connected,  unbalanced  flywheels.  This  thesis  also  provides 
a  flexible  dynamics  model  of  vibrations  that  can  be  used  to  study  various  Integrated 
Power  and  Attitude  Control  System  (IPACS)  configurations. 

1.1  Definitions 

For  the  purpose  of  this  thesis,  a  battery  will  refer  to  a  secondary  electrochemical 
cell  battery.  A  flywheel  is  a  rotating  mass  which  is  used  to  store  kinetic  energy.  An 
advaneed  flywheel  will  be  a  flywheel  unit  consisting  of  a  high-speed,  high-moment  of 
inertia  (MOI),  low- mass  rotor,  frictionless  electromagnetic  bearings,  a  brushless  elec¬ 
tric  motor/generator,  and  the  electronics  necessary  to  control  the  motor  and  bear¬ 
ings,  as  well  as  the  associated  support  structure.  In  this  thesis,  advanced  flywheels 
are  assumed  but  not  explicitly  stated  each  time.  Flywheels  as  studied  here  are  those 
primarily  intended  for  energy  storage,  which  excludes  similar  reaction  wheel  and  con¬ 
trol  moment  gyro  systems.  A  satellite  flywheel  energy  storage  system,  usually  referred 
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to  in  the  context  of  an  IPACS,  contains  a  minimum  of  two  flywheel  units  in  order 
to  allow  for  a  net  zero  angular  momentum  and  prevent  uncontrolled  spinning  of  the 
satellite.  Four  flywheel  units  are  required  for  full,  uncoupled,  3-axis  attitude  control 
and  energy  storage  in  a  non-gimbaled  conhguration. 

Whirl  is  a  natural,  resonant,  rigid  body,  gyroscopic  vibration  mode  in  the  form 
of  a  precession  motion  that  occurs  in  a  rotor/bearing  system.  It  is  described  in  Sec¬ 
tion  2.3.5.  Extra- synchronous  whirl  excitation  refers  to  two  whirl  modes  at  other-than- 
wheel  speeds:  sub- synchronous,  which  is  below  the  speed  of  the  rotor  in  question,  and 
super-synchronous,  which  is  above  it.  A  spinning  unbalanced  rotor  causes  a  vibration 
input  at  its  own  wheel  speed  (spin  speed),  but  there  are  resonant  vibration  modes 
at  frequencies  other  than  the  spin  frequency.  In  a  multiple  flywheel  system,  a  second 
rotor  provides  a  direct  source  of  vibration  at  extra-synchronous  speeds. 

1 . 2  Overview 

All  satellites  have  electronic  equipment  which  requires  electrical  energy  to  run. 
Since  power  from  solar  cells  is  not  available  continuously,  satellites  need  an  energy 
storage  subsystem.  During  periods  of  excess  power  generation,  extra  energy  is  stored 
in  a  battery.  When  the  satellite  needs  more  power  than  the  solar  cells  can  generate,  it 
uses  the  energy  stored  in  the  battery.  One  alternative  to  chemical  batteries  for  energy 
storage  is  advanced  flywheels.  Flywheels  have  not  yet  been  used  for  energy  storage  in 
any  space  missions,  but  the  technology  is  maturing  quickly,  and  someday  they  may 
be  a  viable  alternative  for  the  satellite  designer. 

One  additional  beneht  of  flywheel-based  energy  storage  is  its  inherent  ability  to 
control  the  attitude  of  a  satellite.  Many  satellites  use  some  form  of  momentum  ex¬ 
change  device  for  attitude  control.  Since  a  flywheel  system  has  multiple  rotating  wheels 
it  can  change  the  satellite’s  attitude  by  exchanging  momentum  between  flywheels  and 
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the  spacecraft.  Thus  an  IPACS,  if  well  designed,  can  save  weight  by  combining  two 
necessary  components  of  the  satellite  bus. 

There  are  two  key  performance  measures  for  a  satellite  energy  storage  system: 
high  specihc  energy  and  high  specihc  power.  In  addition,  the  energy  storage  system 
must  meet  several  requirements  in  order  to  be  considered  for  use  in  space.  It  must 
operate  maintenance-free  in  widely-varying  temperature  conditions  and  survive  in  a 
radiation  environment.  It  must  perform  under  these  conditions  throughout  a  long 
lifetime — often  greater  than  10  years,  with  multiple  daily  charge/discharge  cycles 
during  its  entire  operational  lifetime.  It  must  survive  a  harsh  vibration  regime  during 
launch,  without  itself  creating  unwanted  vibration  in  the  satellite  during  operation. 
Finally,  it  must  be  completely  reliable  from  the  beginning  of  the  satellite’s  lifetime  to 
the  end. 

Batteries,  used  for  energy  storage  on  every  satellite,  are  far  from  ideal.  They  can 
provide  either  high  specihc  power  or  high  specihc  energy,  but  not  both.  They  have  a 
limited  lifetime,  measured  not  just  from  the  beginning  of  their  service  life,  but  from 
the  date  of  manufacture.  They  also  require  a  carefully  controlled  thermal  environment 
to  avoid  performance  loss  or  even  damage.  They  do  excel  in  a  few  key  areas,  however. 
They  are  relatively  simple  and  create  zero  external  disturbances.  Most  importantly, 
the  technology  is  mature  and  there  is  a  long  history  of  battery  usage  on  satellites. 
They  have  proven  to  be  predictable  and  reliable  when  used  in  a  well-designed  system. 

Flywheel  energy  storage  systems  have  not  yet  been  used  in  any  space  system, 
but  in  some  ways  their  theoretical  performance  is  far  better  than  that  of  batteries.  A 
hywheel  system  is  able  to  satisfy  demands  for  high  specihc  energy  and  power.  It  can 
theoretically  do  so  without  signihcant  degradation  for  an  extremely  long  life  measured 
both  in  time  and  in  charge/discharge  cycles — lifetime  is  a  minor  design  factor,  but 
some  current  plans  call  for  hywheels  designed  to  operate  for  15  years  and  90  thousand 
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cycles.  A  flywheel  system  can  operate  in  any  thermal  environment  suitable  for  the 
satellite. 

Unfortunately  for  the  satellite  designer,  however,  flywheel  technology  is  far  less 
mature  than  battery  technology.  Flywheel  systems  have  not  yet  been  developed  at 
scales  appropriate  for  satellites.  Current  systems  exist  only  as  bench-top  research 
units,  and  the  power  supplies  and  drive  electronics  have  not  been  scaled  to  an  ap¬ 
propriately  small  size.  Also,  flywheel  systems  have  not  yet  demonstrated  the  required 
reliability.  Furthermore,  rotating  unbalanced  rotors  inherently  create  vibration  which 
must  be  eliminated  or  at  least  mitigated  to  avoid  affecting  satellite  operations  nega¬ 
tively. 

For  attitude  control,  no  direct  comparison  between  batteries  and  flywheels  can 
be  made.  Instead,  flywheels  can  be  compared  to  the  entire  battery  and  momentum 
exchange  attitude  control  systems.  Batteries  are  intended  as  energy  storage  devices 
only,  and  with  no  moving  parts,  they  offer  no  attitude  control.  On  the  other  hand, 
any  advanced  flywheel  system  is  more  than  adequate  to  offer  attitude  control  for  a 
satellite  in  at  least  one  dimension.  Solving  all  of  the  other  problems  of  creating  a 
space-worthy  high  performance  flywheel  will  ensure  that  the  system  is  able  to  control 
the  rotor  momenta  sufficiently  to  orient  the  spacecraft.  Some  minor  concerns  are  a 
slight  oversizing  of  the  system — to  ensure  sufficient  margin  for  both  energy  storage 
and  attitude  control — and  appropriate  geometry  and  control  laws.  The  control  laws 
for  flywheel  attitude  control  are  non-trivial,  but  engineers  have  been  developing  them 
for  decades  and  they  currently  await  hardware  implementation. 

While  there  are  still  hurdles  in  the  way  of  widespread  flywheel  system  adop¬ 
tion  on  satellites,  an  incredible  array  of  challenges  have  already  been  solved.  High 
tensile  strength  carbon  hbers  enable  the  creation  of  light  and  strong  rotors.  Actively 
controlled  contact-free  magnetic  bearings  waste  no  energy  as  heat  due  to  friction. 
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with  minor  and  controllable  system  losses  in  other  areas.  Advanced  brushless  motor/ 
generators  allow  for  very  efficient  energy  storage.  Good  design  and  robust  computer 
control  of  the  bearings  and  motors  enables  stability  throughout  the  operating  regime. 
Two  problems  that  remain  are  scaling  the  technology  to  a  reasonable  size  and  ensuring 
an  acceptable  level  of  system  vibration. 

After  solving  all  technical  problems  satisfactorily,  there  are  two  remaining  steps 
to  be  completed  before  widespread  adoption  of  advanced  flywheels  is  possible.  First, 
a  reasonably  sized,  operationally  representative  unit  needs  to  be  built  and  tested.  To 
date,  most  work  has  been  performed  on  either  larger  systems  or  component-wise  on 
smaller  parts.  A  reasonably  sized  unit  would  include  all  necessary  components  in  a 
package  small  enough  to  £t  in  a  simple  technology  demonstrator.  Finally,  a  successful 
technology  demonstration  satellite  must  be  flown  to  give  other  satellite  designers  proof 
that  flywheels  are  viable  in  space. 

1.3  Objectives 

This  thesis  will  examine  the  problem  of  flywheel-induced  vibrations  on  satellites, 
focusing  on  the  interactions  between  multiple  imperfect  flywheels  at  varying  speeds 
and  the  vibration  inputs  this  imbalance  creates  for  a  satellite.  Even  the  most  precisely 
manufactured  flywheels  have  some  residual  imbalance.  At  wheel  speeds  of  high  tens  of 
thousands  of  revolutions/minute  (RPM),  this  vibration  creates  a  potentially  harmful 
amount  of  vibration.  Bearings  and  soft  mounts  reduce  this  vibration,  but  they  cannot 
completely  eliminate  it. 

In  addition  to  the  individual  flywheels,  the  entire  IPACS  consisting  of  multiple 
flywheels  must  keep  vibration  within  an  acceptable  limit.  The  envelope  must  include 
the  interaction  of  multiple  wheels  operating  at  different  combinations  of  speeds.  Pre¬ 
vious  flywheel  vibration  research  has  focused  primarily  in  single-rotor  vibration  rather 
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than  vibration  interaction.  This  thesis  will  develop  a  linear  state-space  model  to  inves¬ 
tigate  potential  sonrces  of  low  freqnency  system  excitation  caused  by  beat  phenomena 
and  extra-synchronous  whirl  excitation  as  two  connected  flywheels  operate  as  part  of 
an  IPACS.  A  linear  model  is  sufficient  to  prove  the  existence  of  gyroscopic  vibra¬ 
tion  interaction.  The  model  will  be  used  to  study  a  two-flywheel  system.  However, 
the  model  is  flexible  enough  to  be  used  in  future  investigations  of  IPACS  with  an 
arbitrary  number  of  individually  oriented  flywheels. 

This  thesis  will  seek  to  answer  the  question  of  whether  the  beat  frequencies 
caused  by  similar  flywheel  rotation  speeds  or  the  extra-synchronous  interactions  be¬ 
tween  multiple  connected  flywheels  can  cause  harmful  low  frequency  vibration  in  an 
IPACS.  Results  will  be  limited  by  the  model  assumptions:  hxed-satellite,  small  an¬ 
gle  rotations,  linear  springs,  and  limited  geometry  and  input  conhgurations.  These 
assumptions  are  described  in  Section  2.5. 

1.4  Organization 

This  thesis  is  organized  as  follows:  Chapter  I  provides  a  brief  overview  of  the 
thesis.  Chapter  II  reviews  relevant  literature  on  the  subject  of  flywheels  and  provides 
background  information  in  support  of  the  thesis.  The  modeling  methodology  and 
validation  are  covered  in  Chapter  III.  Chapter  IV  details  the  results  of  the  analysis, 
and  Chapter  V  summarizes  the  conclusions  and  recommendations. 
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II.  Background 


2. 1  Literature  Review 

A  flywheel  is  a  device  that  stores  rotational  kinetic  energy  for  later  use  in  the 
form  of  a  rotating  mass.  Flywheels  have  countless  applications,  both  realized  and  the¬ 
oretical,  but  one  as  yet  unrealized  application  is  the  use  of  flywheels  onboard  a  satellite 
for  energy  storage.  A  flywheel  system — at  least  two  wheels  would  be  necessary — could 
supplement  or  even  replace  the  secondary  cell  electrochemical  batteries  on  the  satel¬ 
lite.  In  addition,  with  appropriate  control  algorithms,  the  flywheel  system  could  be 
used  to  provide  both  energy  storage  and  attitude  control  to  the  satellite.  IPACS  offers 
potential  performance  benehts  and  weight  savings  (and  consequently  cost  savings)  to 
the  satellite  designer.  Flywheels  have  practical  potential  applications  on  the  ground 
as  well  as  in  space,  but  this  review  will  be  primarily  limited  to  space  applications. 

As  described  by  Genta,  flywheels  have  been  used  for  millenia,  from  the  inven¬ 
tion  of  spindles  and  potter’s  wheels.  The  high  inertia  of  a  rotating  flywheel  smooths 
out  changes  in  the  motion  of  a  rotating  body.  Motors  often  rely  on  this  smooth¬ 
ing  for  steady  output,  and  many  types  of  motor  would  not  operate  at  all  without  a 
flywheel  (5:3,16). 

Flywheels  have  a  long  engineering  heritage,  but  Sputnik  was  only  launched  in 
1957,  so  the  problem  of  energy  storage  in  space  is  just  over  50  years  old.  A  fly¬ 
wheel  energy  storage  system  for  satellites  was  hrst  proposed  by  Roes  only  a  few  years 
later  (7:17-18). 

No  IPACS  paper  would  be  complete  without  a  reference  to  Roes.  In  1961,  he 
proposed  a  flywheel  system  for  satellite  energy  storage  (7:8).  He  did  not  consider 
using  flywheels  for  attitude  control,  only  for  energy  storage.  The  idea  of  combining 
the  attitude  control  and  energy  storage  functions  of  a  flywheel  into  an  IPACS  began  to 
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appear  in  the  early  1970’s  in  several  technical  papers  from  the  National  Aeronautics 
and  Space  Administration  (NASA)  Langley  Resarch  Center  (7:17-18). 

Early  flywheel  studies  were  not  limited  to  the  investigation  of  advanced  fly¬ 
wheels.  In  a  1997  report,  Hall  states  that,  “several  studies  in  the  1960’s  and  1970’s 
indicated  that  the  use  of  steel  flywheels  on  mechanical  bearings  would  be  competitive 
with  the  chemical  batteries  of  the  time  (7:5).” 

Since  then,  flywheels  have  been  continually  studied  for  space  applications.  NASA 
commissioned  the  enormously  detailed  Integrated  Power/ Attitude  Control  System 
(IPACS)  Study  in  1974  and  a  similarly  thorough  Advaneed  Integrated  Power  and 
Attitude  Control  System  (IPACS)  Study  in  1985.  Contemporary  theoretical  system 
performance  has  continued  to  climb,  but  not  at  the  rate  anticipated  by  some  de¬ 
signers.  Practical  system  performance  by  necessity  lags  even  further  behind.  In  1976 
NASA  predicted  a  system  energy  density  of  300  W-hr/kg  by  the  year  2000.  By  1992, 
this  prediction  was  lowered  to  100  W-hr/kg  with  a  statement  that  not  all  failure 
modes  or  safety  needs  had  yet  been  identified  for  such  systems  (7:10). 

Progress  has  been  made  in  a  variety  of  areas  since  then,  however.  Among  these 
advances  in  theory  and  application  are  high  efficiency  electric  motors,  magnetic  bear¬ 
ings,  composite  rotors,  and  advanced  attitude  control  algorithms.  In  1985,  Genta 
published  the  seminal  Kinetie  Energy  Storage,  which  was  an  attempt  to  summarize 
flywheel  engineering  efforts  for  all  applications  and  provide  an  up  to  date  review  of 
the  subject  (5:v). 

While  much  advanced  flywheel  technology  has  been  developed  for  aerospace  ap¬ 
plications,  the  benehts  are  beginning  to  spread  to  other  industries.  By  1996  advanced 
flywheels  were  beginning  to  be  considered  for  use  in  several  terrestrial  applications, 
including  uninterruptible  power  supplies  and  hybrid  vehicles  (20).  The  U.S.  Navy  is 
currently  using  flywheels  and  associated  technology  in  development  and  helding  of 


their  Electromagnetic  Aircraft  Launch  System  to  replace  steam  catapults  on  aircraft 
carriers  (3). 

Active  magnetic  bearings  are  a  critical  enabling  technology  for  advanced  fly¬ 
wheels.  In  a  vacuum  they  are  frictionless,  reducing  or  eliminating  many  of  the  prob¬ 
lems  otherwise  associated  with  bearings  for  such  high  speed  devices.  The  controllable 
bearing  stiffness  can  be  low,  isolating  the  inherent  vibrations  caused  by  an  unbal¬ 
anced  rotor.  In  addition,  Liters  and  other  control  algorithms  can  be  used  to  control 
the  bearings  such  that  the  the  vibration  isolation  is  tuned  to  problem  areas  (15:2). 

At  the  turn  of  the  century,  flight  prospects  for  flywheels  looked  great.  NASA  had 
been  working  with  with  U.S.  Flywheel  Systems,  Inc.,  TRW,  Texas  A&M  University, 
the  University  of  Texas,  and  Boeing  for  hve  years  to  develop  and  build  advanced 
flywheels.  These  efforts  were  rewarded  in  December  of  1999  with  a  successful  test  run 
of  their  D1  flywheel  unit  at  60,000  RPM — a  then-world  record  for  a  magnetic  bearing 
flywheel.  These  efforts  were  intended  to  lead  to  a  technology  demonstration  payload 
for  the  International  Space  Station  (ISS). 

By  2001,  plans  were  firmly  in  place  for  the  most  promising  space-based  technol¬ 
ogy  demonstration  to  date.  NASA  had  continued  to  work  on  plans  for  a  technology 
demonstration,  and  onboard  testing  of  a  Flywheel  Energy  Storage  System  for  the  ISS 
was  to  begin  by  2005.  A  successful  test  of  the  system  could  have  led  to  the  eventual 
replacement  of  the  station’s  batteries  (14:2).  Unfortunately,  funding  for  this  program 
was  cancelled  in  2002  (2),  and  no  further  solid  plans  for  a  technology  demonstration 
have  been  made. 

Some  residual  flywheel  development  continued,  however.  In  July  of  2003,  NASA 
demonstrated  a  basic  IPACS  capability  on  an  air  bearing  table  with  two  flywheel 
modules  (19:64).  The  test  setup  can  be  seen  in  Figure  1.  Meanwhile,  NASA  was 
working  on  the  more  advanced  G2  flywheel  module  with  better  performance.  This 
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was  the  first  flywheel  module  designed  in-house  by  NASA’s  Glenn  Research  Center. 
It  had  a  higher  power  rating,  lower  spin  losses,  and  a  more  generous  thermal  envelope 
than  the  previous  D1  flywheel. 


Figure  1:  Photo  of  NASA’s  D1  and  HSS  flywheels  demonstrating  integrated  power 
and  attitude  control,  July  2003  (19:64) 

In  September  of  2004,  NASA  successfully  tested  the  G2  flywheel  module  to  a 
speed  of  41,000  RPM.  G2  is  shown  in  Figure  2  (9:132).  Later  that  month,  the  same 
team  placed  two  flywheel  modules  (the  older  D1  and  the  newer  G2)  on  an  air  bear¬ 
ing  table  to  demonstrate  a  full  IPACS  capability.  They  succeeded  in  demonstrating 
controllable  torque  up  to  ±0.8  N-m  and  power  transfer  from  0-300  W.  This  was  a 
hrst  for  high-power,  high-speed  IPACS.  Since  then,  however,  active  efforts  towards  a 
flight-worthy  technology  demonstration  have  remained  stalled. 

With  more  recent  technology  advances,  some  researchers  are  beginning  to  de¬ 
sign  small  IPACS  for  small  satellites  with  demanding  mission  prohles.  Lappas  et  al. 
discussed  this  in  an  article  that  appears  to  be  the  most  recent  comprehensive  review 
of  IPACS  history  and  literature  (12). 

The  control  algorithms  for  IPACS  have  been  studied  very  thoroughly,  with  many 
papers  written  about  both  energy/power  storage  and  attitude  control.  Two  types  of 
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Figure  2:  Photo  of  NASA’s  G2  flywheel  module,  tested  to  41  kRPM  in  2004  (10:95) 


attitude  control  models  have  been  studied:  gimbaled  CMG-like  flywheels  and  hxed 
momentum  wheel-like  flywheels  (7:22).  Hall  demonstrated  that  with  the  minimal  mo¬ 
mentum  wheel  conhguration  (4  flywheels),  the  attitude  control  and  energy  storage 
functions  of  an  IPAGS  system  could  be  completely  decoupled.  This  simplihes  control 
development  for  both  functions  (6:1894). 

To  date,  the  study  of  IPAGS  vibrations  has  been  very  limited.  Previous  IPAGS 
research  has  focused  heavily  on  control  algorithms,  assuming  rigid  bearings  and  per¬ 
fectly  balanced  flywheels.  In  his  Ph.D.  dissertation.  Park  studied  an  IPAGS  with 
flexible  magnetic  bearings,  unbalanced  flywheels,  and  flexible  appendages.  He  used 
the  imbalances  and  appendages  as  inputs  to  an  IPAGS  to  develop  control  algorithms 
and  physical  means  to  mitigate  problems  in  an  IPAGS.  He  found  that  wheel-speed 
notch  hlters  in  the  control  algorithm  and  vibration  control  masses  on  the  end  of  flex¬ 
ible  satellite  components  could  effectively  reduce  vibration  and  power  surge  problems 
in  the  IPAGS  and  consequently  in  the  satellite  (15).  Park’s  physical  model  was  similar 
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to  the  one  developed  in  this  thesis,  but  this  thesis  investigates  instead  the  interactions 
of  vibrations  between  multiple  flywheels. 

The  beat  frequency  has  been  studied  briefly  in  rotordynamics.  Research  per¬ 
formed  at  the  Naval  Postgraduate  School  used  the  beat  frequency  as  a  tuning  aid  to 
match  hlter  frequencies  to  vibration  frequencies  (13).  The  equations  of  motion  for  a 
rigid  body  gyroscope  are  linear,  however,  and  since  the  beat  frequency’s  effects  on 
a  linear  system  are  small  as  shown  in  Section  2.3.6,  there  has  not  been  much  study 
of  the  effects  of  beat  frequency  in  flywheel  systems.  One  exception  is  a  NASA  in¬ 
vestigation  on  beat  frequency  effects  caused  by  pulse  width  modulation  of  position 
sensor  signals  for  a  magnetically  suspended  flywheel  rotor,  but  this  model  dealt  with 
non-linear  effects  rather  than  linear  physical  behavior  (11). 

While  IPACS  vibration  has  been  largely  neglected,  the  study  of  a  single  fly¬ 
wheel’s  vibration  is  completely  within  the  realm  of  rotordynamics,  which  is  very  well 
studied  and  is  applicable  to  a  wide  variety  of  modern  mechanical  systems.  Vance’s 
Rotordynamics  of  Turbomachinery  is  one  of  the  early  comprehensive  books  on  the  sub¬ 
ject  (22:iv).  Research  in  the  held  of  rotordynamics  is  ongoing,  and  with  the  advent  of 
realizable  controlled  magnetic  bearings,  the  state  of  the  art  continues  to  advance. 

IPACS  research  has  been  ongoing  for  50  years  now,  and  no  immediate  demon¬ 
stration  is  planned.  The  main  obstacle  in  the  way  of  a  successful  hight  demonstration 
is  the  complexity  the  problem — using  an  advanced  and  dynamic  system  to  perform 
two  unrelated  tasks.  While  many  of  these  complications  have  been  studied  and  some 
of  them  have  been  mitigated,  the  vibration  interactions  between  multiple  hywheels  in 
one  IPACS  have  been  neglected.  This  thesis  seeks  to  provide  a  look  into  the  problem 
of  vibration  interaction  in  order  to  determine  whether  it  will  be  a  problem  for  IPACS 
in  space.  To  peform  this  study,  a  hexible  dynamics  model  is  developed  that  will  enable 
research  into  vibrations  of  multiple  gyroscope  systems. 
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2.2  Coordinates  and  Nomenclature 


There  are  two  types  of  coordinates  used  in  this  model:  wheel-aligned  and  global. 
Both  coordinate  frames  are  right-handed  in  x,  y,  and  z,  with  corresponding  rotations 
and  0;^. 

The  wheel-aligned  coordinates  are  oriented  such  that  the  2;  axis  is  along  the 
flywheel  spin  axis  for  the  wheel  in  question.  The  spin  axis  always  points  away  from 
the  support  structure,  which  is  the  system  center  of  mass  (COM).  The  x  and  y  axes 
are  arbitrary  in  this  coordinate  frame  since  all  inputs  and  outputs  are  cyclical  in 
nature  about  z.  All  discussions  of  individual  flywheels  refer  to  these  local  coordinates, 
including  all  imbalance-induced  input  forces  and  wheel  speeds. 

The  global  coordinates  are  arranged  with  z  “up”.  This  coordinate  frame  is 
arbitrarily  aligned.  All  system  level  references — including  all  response  plots  shown  in 
this  thesis — will  use  the  global  coordinate  system.  Model  inputs  are  applied  internally 
in  the  global  frame.  Coordinate  systems  are  arranged  as  shown  in  Figure  3. 


(a)  Global  Coordinates  (b)  Wheel-aligned  Coordinates 

Figure  3:  Coordinate  systems  used  in  this  thesis 


Nomenclature  is  dehned  as  it  is  introduced,  as  well  as  appearing  in  a  nomen¬ 
clature  section  at  the  beginning  of  this  thesis.  A  few  terms  are  potentially  confusing, 
so  a  preview  of  them  is  in  order  here.  Rotation  angles  Ox,y,z  are  defined  as  shown  in 
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Ox  Oy  Oz 


,  which  are  rotations 


Figure  3.  6  is  the  vector  of  rotational  states 
about  X,  y,  and  z,  respectively.  In  the  discussion  of  flywheels  and  rotors,  wheel  rota¬ 
tion  and  speed  are  of  particular  importance,  so  0  with  no  letter  subscript  will  refer  to 
the  rotation  angle  of  a  rotor.  Similarly,  u  with  no  letter  subscript  will  refer  to  rotor 
speed,  which  would  otherwise  be  referred  to  as  0^- 


Furthermore,  there  are  three  sets  of  units  in  widespread  use  to  describe  rota¬ 
tional  speed:  rad/s,  RPM,  and  Hz  or  revolutions/second  (RPS).  Flywheel  dynamics 
must  be  calculated  in  radians,  but  discussion  is  more  intuitive  in  units  of  Hz  or 
RPS.  Finally,  system-level  discussions  of  high-speed  advanced  flywheels  commonly 
uses  units  of  RPM.  The  model  developed  in  this  thesis  uses  units  of  radians  and 
rad/s  internally,  but  most  discussion  of  wheel  speeds  in  this  paper  will  be  in  terms  of 
Hz  and  RPS. 


2. 3  Fundamental  Equations 

Understanding  of  several  basic  sets  of  equations  is  required  for  the  study  of 
flywheel  motion.  A  quick  overview  of  some  of  these  concepts  is  provided  here. 

2.3.1  Equation  of  Motion  for  a  Gyroscopic  Body.  In  this  model,  flywheels  are 
modeled  as  gyroscopic  rigid  bodies.  The  equation  of  motion  (EOM)  for  a  gyroscopic 
rigid  body  is  shown  in  matrix  form  in  Equation  1  (18:124).  Equation  1  also  describes 
the  equation  of  motion  for  a  non-spinning  rigid  body  since  gyroscopic  stiffness,  G,  is 
a  function  of  wheel  speed,  u,  and  it  is  zero  for  a  non-spinning  body. 


Mq(t)  +  {C  +  G  (oj(t)))  q(t)  +  Kq(t)  =  u(t) 


(1) 


where  the  state  vector  q  = 


is  composed  of  r  = 


[X  y  z 


and  6  = 


[dx  Oy  6z]^  1  is  a  vector  of  input  forces,  M  is  the  mass  matrix  described  in  terms 
of  body  mass  m  and  directional  MOI  Ix,y,z-  If  the  body  is  represented  as  a  point 


mass,  M  is  the  diagonal  matrix  M  =  diag 


m  m  m  Ix  ly  Iz 


.  K  and  C  are 


the  stiffness  and  damping  matrices  between  the  body  and  external  nodes,  and,  when 
described  in  wheel-aligned  coordinates,  G  is  the  sparse  skew-symmetric  gyroscopic 
matrix  shown  below. 


G  = 


0  —Iz(^  0 

IzCO  0  0 

0  0  0 


In  the  simple  case  of  a  single  body  connected  by  springs  to  a  hxed  support. 


stiffness  would  be  written  as  K  =  —diag 


kx  ky  kz  Ke,  Ke„ 


with  trans¬ 


lational  and  rotational  stiffnesses  k  and  k,  respectively.  Likewise,  damping  of  a  single 


body  would  be  described  by  C  =  —diag 
lational  and  rotational  damping  of  c  and  C. 


Cx  C^i  Cz  Cg 


with  trans- 


In  this  model,  however,  the  bodies  will  be  connected  not  to  hxed  supports,  but 
to  each  other.  Proper  modeling  of  this  inter-body  stiffness  requires  the  use  of  different 
(and  more  complicated)  stiffness  and  damping  terms.  These  inter-body  stiffness  and 
damping  terms  are  only  applicable  to  a  system  rather  than  to  individual  rigid  bodies, 
and  they  will  be  discussed  in  Section  3.3.1. 


2.3.2  State-Space  Equation  of  Motion.  State-space  representation  is  a  con¬ 
venient  format  for  writing  and  solving  linear  differential  equations,  and  it  is  the  form 
that  will  be  used  for  the  model  in  this  thesis.  In  general,  a  linear  system  can  be 
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described  in  state-space  according  to  Equation  2  (17:210).  Note  that  the  notational 
dependence  on  time  is  dropped  for  convenience. 


q  =  Aq  Bw  (2) 

where  A  is  the  state  matrix  that  describes  system  behavior  and  B  is  an  input  matrix 
linking  input  forces  to  states. 

Equation  3  below  shows  Equation  1,  the  gyroscopic  EOM,  in  state-space  form 
with  the  addition  of  an  input  matrix,  B.  Note  that  the  entire  matrix  EOM  is  found  in 
the  second  half  of  this  equation;  the  top  half  is  simply  a  computational  convenience. 


q 

0  I 

q 

+ 

0 

q 

M-^K  M-i(C  +  G) 

q 

- 1 

PQ 

rH 

1 

2.3.3  Rotated  Equation  of  Motion.  If  the  flywheel  EOM  is  known  in  the 
wheel-aligned  coordinate  frame,  but  the  integration  is  to  be  carried  out  in  a  different 
global  coordinate  frame.  Equation  3  must  be  rotated  accordingly.  Recalling  that  the 
entire  matrix  EOM  is  found  in  the  second  half  of  the  state-space  equation,  the  correct 
application  of  the  rotation  matrix  R  is  shown  below  in  Equation  4.  Proper  rotation  is 
required  when  a  collection  of  bodies  with  different  local  coordinate  frames  is  integrated 
into  one  state-space  system. 


q 

I 

0 

q 

0 

R 

0  I 

M-^K  M-i(C  +  G) 


I 

0 

q 

+ 

0 

RT 

q 

0 

RM-^BRT 


u  (4) 
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2.3.4  Centripetal  Force.  The  primary  input  in  this  model  will  be  the  vibra¬ 
tion  caused  by  an  unbalanced  spinning  flywheel  rotor,  which  is  a  centripetal  force. 
Equations  5  and  6  show  the  equations  for  centripetal  acceleration,  ttc,  and  its  associ¬ 
ated  force,  fc  in  polar  coordinates  (1:75-76). 


ttc  =  —pu^ 

(5) 

m2  2 

jc  =  —euj  =  —mpoo 

(6) 

where  eccentricity  e  represents  rotor  mass  m  and  the  distance  between  the  flywheel’s 
COM  and  the  center  of  the  shaft,  p.  Wheel  speed  is  again  represented  by  u.  The 
distance  p  is  shown  in  Figure  4.  The  vectors  e  and  p  point  from  the  center  of  mass 
to  the  flywheel  shaft. 


Figure  4:  A  spinning  flywheel  of  mass  m  with  COM-shaft  distance  p  creates  a  cen¬ 
tripetal  input  force  proportional  to  eccentricity  e  =  mp 


2.3.5  Natural  Frequencies  of  a  Rotor/Bearing  System.  The  IPACS  model  in 
this  thesis  will  be  used  to  evaluate  vibration  interactions.  One  potentially  troublesome 
vibration  is  resonant  frequency  excitation  of  the  flywheels.  Rigid  rotors  have  several 
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vibration  modes,  and  those  are  briefly  explained  here.  In  this  brief  discussion  of  natural 
frequencies,  a;i-a;4  will  refer  to  natural  frequencies  of  the  rotor/bearing  system,  and 
u  with  no  subscripts  will  refer  to  wheel  speed. 

When  a  rigid  rotor  is  restrained  in  the  x  and  y  directions  by  springs  of  stiffness 
K  as  shown  in  Figure  5,  the  EOM  from  Equation  1  can  be  simplified  to  the  forms 
shown  below  in  Equations  7-10. 


mx  +  2kx  =  0  (7) 

my  +  2ky  =  0  (8) 

It^x  +  IpX)9y  +  -kL‘^9x  =  0  (9) 

lT9y  —  Ipxj9x  +  =  0  (10) 


where  m  =  rotor  mass,  L  =  rotor  length,  and  It  and  Ip  are  the  transverse  and  polar 
MOIs,  respectively  (22:125). 

The  natural  frequencies,  Un,  of  the  rigid  body  rotor/bearing  system  shown  in 
Figure  5  and  described  in  Equations  7-10  are  shown  in  Equations  11-13.  Recall  that 
u  represents  angular  rotor  speed. 


Ui  =  U2  = 


us{u) 

U4^{u) 


2K 
m 

Ip 


_  Ip 


IKL^ 


2h 


+ 


'KL^ 

2It 


2h 


-u 


2h 


-u 


(11) 

(12) 

(13) 


Natural  frequencies  Ui  and  1x2  correspond  to  translational  vibration  modes  in 
the  X  and  y  directions.  Erequencies  cus  and  0)4  correspond  to  the  forward  and  backward 
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Figure  5:  A  long  rigid  rotor  constrained  by  springs  in  x  and  y 


whirling  modes,  respectively.  Whirling  is  a  precession-like  motion.  Forward  whirl  (cua) 
is  rotation  in  the  same  direction  as  wheel  spin,  as  shown  in  Figure  6.  Backward  whirl 
{(jJi)  is  rotation  in  the  opposite  direction. 

Modes  1  and  2  are  constant  and  identical.  They  represent  “bouncing”  in  the 
X  and  y  directions.  There  is  no  angular  motion  associated  with  these  modes.  The 
frequency  of  these  modes  is  solely  a  function  of  bearing  stiffness  and  rotor  mass. 


Modes  3  and  4  are,  respectively,  the  speed-dependent  forward  and  backward 
conical  whirling  motions.  With  a  wheel  speed  of  zero,  a;3(0)  =  a;4(0)  =  ut  =  \J 
which  is  the  natural  rigid  body  pitching/yawing  frequency.  Rigid  body  vibration  at 
this  frequency  will  be  simple  pitching  or  yawing  about  9^  or  6y. 


As  rotor  speed,  ca,  increases,  the  pitching  and  yawing  motions  become  whirling 
motions,  and  the  frequencies  of  the  forward  and  backward  whirls  diverge  from  cot,  with 
forward  whirl  speed  increasing  and  backward  whirl  speed  decreaseing.  As  wheel  speed 
approaches  inhnity,  approaches  P  and  a;4  approaches  0.  P  is  the  ratio  Ip/ It  (polar 
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Figure  6:  Depiction  of  forward  whirling  motion 

MOI/transverse  MOI),  which  represents  rotor  conhguration.  For  a  short  coin-like  disc 
or  hoop,  P  =  2.  For  an  inhnitely  long  rotor,  P  =  0  (22:126). 

Figure  7  shows  the  normalized  natural  frequencies  of  the  forward  and  backward 
whirl  modes,  respectively.  Recall  that  ut,  the  natural  rigid  body  pitching/yawing 
frequency  is  ^ Backward  whirl  has  a  negative  frequency  value  because  it  is 
in  the  direction  opposite  wheel  spin.  These  modes  are  dependent  on  P,  the  rotor 
conhguration  (22:127). 

Of  the  whirl  modes,  cus  is  the  mode  most  easily  excited  by  rotor  imbalance. 
Figure  7(a)  shows  that  short  rotors  are  not  self-exciting  for  forward  whirl,  but  long 
rotors  have  a  critical  coning  frequency  where  the  wheel  frequency  is  synchronous  with 
the  whirl  frequency.  The  critical  coning  frequency,  Ucon  is  found  by  setting  a;  =  cus  in 
Equation  12,  resulting  in  Equation  14. 


^con  — 


(14) 
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(a)  Forward  Whirl  {0J3) 


(b)  Backward  Whirl  (^4) 

Figure  7:  Normalized  whirl  modes  for  various  flywheel  rotor  conhgurations.  Dotted 
line  for  forward  whirl  shows  where  whirl  speed  is  synchronous  with  wheel 
speed.  Natural  rigid-body  pitching  frequency  ut  is  used  to  normalize  plots 
for  all  rotor/bearing  conhgurations  (22:128) 
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Operation  near  the  critical  coning  frequency  can  be  unstable,  although  momen¬ 
tarily  passing  through  this  region  is  acceptable,  allowing  operation  in  supercritical 
regions  (21:354).  For  long  rotors  {P  <  1),  this  critical  frequency  must  be  kept  low 
enough  to  be  out  of  the  operating  range  of  the  rotor.  Shorter  rotors  {P  >  1)  do  not 
have  this  problem,  but  an  outside  excitation  at  the  coning  frequency  could  be  haz¬ 
ardous.  Outside  excitation  at  critical  coning  frequencies  is  examined  in  this  paper  for 
both  short  and  long  rotors,  which  are  the  sub-  and  super-synchronous  rotor  frequency 
problems. 

Natural  frequencies  Un  as  found  in  Equations  11-13  are  described  in  units  of 
rad/s.  For  the  remainder  of  this  thesis,  natural  frequencies  will  be  instead  described 
in  units  of  Hz  as  /„,  where 


_  ^ 
“  27r 


(15) 


This  terminology  will  eliminate  any  confusion  between  wheel  speeds  (described  as  u) 
and  the  natural  frequencies. 


2.3.6  Beat  Frequency.  For  linear  systems,  the  response  of  a  system  to  two 
inputs  can  be  found  by  adding  the  system’s  response  to  the  individual  inputs.  This 
superposition  can  yield  a  stronger  or  weaker  response  than  either  signal  individually, 
depending  on  whether  the  sum  of  the  responses  is  positive  or  negative,  creating  either 
constructive  or  destructive  interference. 

When  waveforms  of  two  frequencies  (cui,  UJ2)  differ  in  frequency  by  a  small 
amount  <5,  {uj2  =  cui  -|-  5),  the  wave  which  results  from  their  combination  will  exhibit 
what  is  known  as  a  beat  phenomenon  due  to  alternating  constructive  and  destructive 
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interference.  Using  the  relationship 


.  „  A+B  A-V 

sin  A  +  sin  B  =  2  sin - cos - , 

2  2 


the  sum  {yi  +  1/2)  of  two  displacements 


yi  =  sincuif 

y2  =  sina;2t  =  sin  (cui  +  S)t 
can  be  written  as  shown  in  Equation  16  (16:23). 


y  =  yi  +  y2 


2  sin 


(16) 


Figure  8  shows  an  example  of  a  beat  frequency  created  by  the  superposition  of 
two  waves  with  similar  frequencies.  The  resultant  waveform  is  a  sine  wave  of  frequency 
u  =  which  is  shaped  by  the  envelope  described  by  ±2  cos  (f)t.  This  envelope 

has  a  much  lower  frequency  than  either  of  the  original  waves. 

When  two  signals  are  combined  in  a  linear  system  they  are  superimposed  addi- 
tively.  In  order  to  see  system  excitation  at  the  beat  frequency,  there  would  have  to  be 
some  non-linear  process  allowing  the  input  frequencies  to  be  multiplied.  Otherwise, 
there  will  be  no  change  in  frequency.  This  lack  of  excitation  at  the  beat  frequency 
should  manifest  itself  in  a  lack  of  energy  at  the  beat  frequency  in  a  power  spectral 
density  (PSD)  plot.  The  signal  shown  in  Figure  8  was  analyzed  with  the  Matlab® 
pwelch  command  to  create  the  PSD  plot  shown  in  Figure  9. 

As  expected.  Figure  9  shows  that  the  signal  in  Figure  8  has  energy  at  both  input 
frequencies,  but  there  is  no  system  energy  at  the  beat  frequency  of  5  Hz.  In  linear 
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Figure  8:  Beat  phenomenon  created  by  inputs  of  frequency  30  and  35  Hz 


Figure  9:  Power  spectral  density  plot  of  the  signal  shown  in  Figure  8.  Note  the  lack 
of  energy  at  the  beat  frequency  of  5  Hz 
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systems  there  should  be  no  excitation  at  this  beat  frequency,  even  if  some  connected 
structure  has  a  natural  frequency  close  to  the  beat  frequency. 

2.Jf.  Model 

The  model  used  in  this  thesis  represents  two  advanced  flywheels  mounted  to 
a  support  structure  and  nominally  spinning  in  opposite  directions.  The  rotors  are 
mounted  with  active  magnetic  bearings,  and  their  operating  speed  range  is  20,000- 
60,000  RPM  (333-1000  RPS).  They  are  powered  by  high-efficiency  motor/generators. 
The  rotors  nominally  spin  at  the  same  speed  to  store  energy,  changing  speed  relative 
to  each  other  to  control  the  attitude  of  the  satellite. 

2. 5  Scope 

The  model  developed  in  this  thesis  relies  on  several  assumptions  to  limit  the 
scope  of  the  analysis.  The  only  source  of  vibration  studied  is  rotor  imbalance.  There 
are  multiple  other  real  sources  of  vibration  including  torque  ripple,  sensor  error,  band¬ 
width  limitations,  and  external  vibrations.  These  additional  vibration  sources  are  ig¬ 
nored. 

Also,  this  model  does  not  take  into  consideration  any  motion  of  the  satellite. 
In  a  satellite  with  an  IPACS,  the  satellite  body  will  be  free  to  rotate,  and  can  be 
controlled  by  the  varying  rotation  rates  of  the  flywheels.  In  the  model  used  in  this 
thesis,  the  IPACS  is  subjected  mainly  to  symmetric  or  periodic  disturbing  forces  from 
a  static  equilibrium  state,  and  it  does  not  experience  large  rotations.  When  necessary, 
a  spring  is  used  to  enforce  small  angles.  The  spring  constraint  makes  the  use  of  small 
angle  approximations  for  IPACS  rotation  appropriate  and  allows  for  a  simpler  linear 
analysis.  Similarly,  this  model  does  not  acconnt  for  any  rotor  translation  along  or 
rotation  about  the  axis  of  the  flywheel,  except  that  the  gyroscopic  stiffness  increases 
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with  increasing  rotor  speed.  The  bearings  in  this  model  are  assumed  to  be  linear 
springs,  as  opposed  to  the  controllable  magnetic  bearings  of  an  actual  IPACS. 

Even  the  best  control  model  will  be  unable  to  completely  hlter  out  all  dis¬ 
turbances  such  as  rotor  imbalance  due  to  limitations  such  as  signal  bandwidth  and 
control  saturation.  This  model  assumes  a  small  residual  amount  of  rotor  imbalance 
that  cannot  be  hltered  out  and  examines  the  interactions  between  multiple  residual 
imbalances.  Therefore,  the  input  and  output  forces  are  small. 

The  scope  of  this  analysis  is  also  limited  to  rigid  flywheel  rotors.  There  are 
higher  bending  modes  associated  with  flexible  rotors,  but  the  first  four  vibration 
modes  discussed  in  Section  2.3.5  are  dominant — bending  modes  are  typically  above  a 
frequency  of  1  kHz  (15:35). 
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III.  Methodology 


3. 1  Overview 

This  thesis  will  use  an  analytical  model  to  study  the  vibration  interaction  of 
multiple  gyroscopes.  Vibrations  can  come  from  several  sources,  but  this  thesis  will 
examine  only  those  caused  by  unbalanced  flywheels.  The  model  will  be  numerically 
integrated  using  the  ode45  command  in  Matlab®. 

3.2  Model  Description 

3.2.1  Model  Construction.  The  model  used  in  this  thesis  is  a  system  of  two 
flywheels  and  a  support  structure  as  shown  in  Figure  10.  The  flywheels  are  axially 
aligned,  with  opposite  spin  directions.  The  flywheels  are  arranged  such  that  the  system 
COM  and  the  support  COM  are  co-located. 


Figure  10:  Basic  conhguration  of  the  model  used  in  this  thesis 

The  flywheels  are  each  connected  to  the  support  structure  with  two  magnetic 
bearings,  as  shown  in  Figure  11.  Flywheels  of  length  I  and  radius  r  are  located  at 
distance  d  from  the  support  COM  and  supported  at  each  end  by  a  magnetic  bearing. 
The  magnetic  bearings  have  only  translational  stiffness,  but  having  one  of  them  at 
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each  end  of  the  rotor  will  create  an  effective  rotational  stiffness.  The  stiffnesses  in 
the  transverse  {x  and  y)  and  axial  {z)  directions  are  separate  and  not  necessarily 
related.  In  this  thesis,  however,  axial  displacements  are  ignored,  so  the  axial  stiffness 
is  unimportant.  When  damping  is  accounted  for  in  this  thesis,  all  springs  shown 
represent  bearings  with  both  stiffness  and  damping. 


r 


Figure  11;  Flywheel  in  housing  connected  to  IPACS  support  structure 

The  model  is  simplihed  by  replacing  each  body  with  a  point  mass  as  shown 
in  Figure  12.  There  is  one  4  degree  of  freedom  (DOF)  spring  (and  damper)  located 
at  the  flywheel’s  COM,  which  is  attached  to  the  support  structure  with  a  rigid  link. 
The  spring  shown  has  transverse  translational  {x  and  y)  and  rotational  {6x  and  6y) 
stiffness. 

For  a  flywheel  of  length  /  with  individual  magnetic  bearing  stiffness  kmag,  the 
model  will  have  linear  stiffness  kmodei  =  ‘^^kmag  and  transverse  rotational  stiffness 
^T, model  ^kmagl-  (22.125).  Damping  is  similar.  Cmodel  ‘^Cmag  and  2^magl  ■ 


j :  wheel  A 


z.O^ 

A 


/:  support  ^ 

Figure  12:  Model  of  support/flywheel  connection.  Bodies  are  point  masses  separated 
by  distance  d 

The  two  flywheels  are  attached  to  the  support  structure  as  discussed  previously 
and  illustrated  again  in  Figure  13.  Also  shown  is  a  connection  to  the  satellite,  which 
represents  a  soft  mount  between  the  support  structure  and  the  rest  of  the  satellite. 
Since  all  forces  studied  are  periodic  or  symmetric  about  the  system  and  support  struc¬ 
ture  COM,  the  support  will  primarily  experience  rotations  rather  than  translations. 
For  this  reason,  the  satellite/support  spring  is  modeled  as  a  3  DOF  spring  with  only 
rotational  stiffness.  The  satellite  in  this  model  is  assumed  to  be  heavy  enough  that 
it  can  be  considered  fixed.  For  most  validation  runs,  the  satellite/support  spring  was 
turned  off  to  allow  free  rotation  of  the  satellite. 

Finally,  an  appendage  can  be  added  to  the  model.  The  appendage  represents  a 
flexible  spacecraft  structure  such  as  a  solar  array  or  antenna,  and  it  is  used  to  study 
low  frequency  excitation.  The  appendage  is  shown  in  Figure  14.  The  appendage  and 
support  structure  COMs  are  co-located,  and  they  are  connected  by  a  2  DOF  spring 
with  only  transverse  {6x  and  6y)  rotational  stiffness. 

3.2.2  Model  Inputs.  The  sources  of  vibration  in  this  model  will  be  rotor 
imbalances.  Real  rotors  can  have  very  small  imbalances  if  they  are  manufactured 
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Figure  13:  System  model  with  satellite  included.  The  satellite/support  spring  can  be 
turned  off  if  needed 


Figure  14:  Model  of  appendage,  which  is  connected  to  the  support  structure  with  only 
a  rotational  spring 
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with  tight  tolerances,  but  imbalances  will  always  be  present  after  manufacturing.  The 
rotors  in  this  model  will  be  assumed  to  be  imbalanced  in  such  a  way  that  they  cause 
a  purely  two-dimensional  vibration:  linear  (in  x,y)  or  rotational  (in  dx,0y).  These 
imbalances  are  shown  conceptually  in  Figure  15.  The  axially  symmetric  imbalances 
in  a  real  rigid  rotor  can  be  described  as  a  combination  of  these  two  imbalances,  but 
they  will  be  examined  individually  in  this  model. 


(a)  Linear  (b)  Rotational 

Figure  15:  Two  sources  of  axially-symmetric  imbalance-induced  vibration 

The  rotating  imbalance  creates  a  centripetal  force  as  discussed  in  Section  2.3.4. 
For  this  model  the  rotor  eccentricities  are  replaced  with  an  ideal  rotor  plus  periodic 
input  forces  synchronized  with  wheel  position  and  proportional  to  eccentricity  and 
the  square  of  the  wheel  speed. 

3.2.3  Model  Parameters.  The  flywheel  for  this  model  is  a  theoretical  flywheel 
only,  but  it  is  intended  to  be  realistically  sized.  Flywheel  parameters  are  given  in 
Table  1. 

Magnetic  bearing  stiffness  and  damping  are  similar  to  those  used  in  some  NASA 
studies  (4).  The  nominal  mass  properties  are  described  in  Hibbeler’s  text  (8).  They 
are  dehned  here  in  such  a  way  that  they  describe  a  short  rotor  (^P  =  >  1  j  to  study 

super-synchronous  whirl.  When  necessary.  It  is  changed  to  adjust  the  rotor  parameter 
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Table  1:  Flywheel  model  parameters 


mag  bearing  stiffness 
mag  bearing  damping 
rotor  mass 
rotor  length 
support/rotor  distance 
rotor  radius 

rotor  shaft/COM  distance 
model  translational  stiffness 
model  rotational  stiffness 
model  translational  damping 
model  rotational  damping 
transverse  MOI 
polar  MOI 
rotor  MOI  ratio 


k 

'^mag 

1756 

kN/m 

Cmag 

3.512 

kN/m/s 

m 

10 

kg 

1 

20 

cm 

d 

15 

cm 

r 

15 

cm 

P 

0.01 

nm 

kmodel 

3512 

kN/m 

^model 

35.12 

kN-m/rad 

^model 

7.024 

kN/m/s 

Cffiodel 

70.24 

N-m/rad/s 

It 

0.0896 

kg-m^ 

Ip 

0.1125 

kg-m^ 

P 

1.2558 

- 

P.  This  adjustment  is  made  to  study  longer  rotors  when  looking  for  sub-synchronous 
whirl.  The  COM-shaft  offset  distance,  p  was  chosen  to  give  similar  disturbance  inputs 
to  the  residual  disturbances  found  by  Park  (15:87).  The  mass  properties  of  the  support 
structure  are  shown  in  Table  2. 

Table  2:  Support  structure  parameters 

mass  m  10  kg 

MOI  10  kg-m^ 


3.3  System  Equation  of  Motion 

3.3.1  System  Equation  of  Motion  Components.  When  a  system  of  indepen¬ 
dent,  unconnected,  gyroscopic  rigid  bodies  is  described  in  state-space  as  shown  in 
Equation  2,  it  takes  the  form  of  the  block  diagonal  matrix  shown  in  Equation  17. 
The  stiffness  and  damping  terms  here  represent  each  body  being  connected  to  a  fixed 
body. 
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where  the  state  vector  q  =  ^  ^  ^  ^  q.  =  ^.t  ^.t  ^,t  ^.t  ^  and  Aj 

and  Bj  are  defined  according  to  Equations  18  and  19,  similar  to  Equation  4.  r*  and 
9i  are  the  position  (x,  y,  z)  and  rotation  {Ox,  Oy,  0^)  vectors  of  the  bodies,  respectively, 
and  Vi  and  uJi  are  the  corresponding  translational  and  rotational  velocities. 


A, 

B. 


0 

m-^k, 

0 

Mr^B, 


I 

Mri(Q  +  G,) 


(18) 

(19) 


A  system  of  connected  rigid  bodies  is  described  by  the  same  equation  of  motion 
given  in  Equation  17  plus  the  addition  of  non-block  diagonal  stiffness  terms  in  the 
A  matrix.  In  the  case  of  this  model,  the  stiffness  between  each  flywheel  and  the 
support  structure  is  identical  in  the  local  axially  aligned  frame.  This  common  stiffness 
matrix  can  be  derived  using  the  spring  equation  with  the  help  of  a  diagram  of  system 
displacements.  With  the  support  and  the  wheel  modeled  as  point  masses,  a  simple 
diagram  of  displacements  in  x  is  shown  in  Figure  16,  which  is  a  simplification  of  the 
housing  support  structure.  Small  angles  are  assumed. 

Figure  16  is  used  to  determine  the  inter-body  spring  forces  caused  by  system 
displacements  according  to  F  =  —kx.  The  bodies  are  modeled  as  point  masses.  Angles 
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j :  wheel 

z,e^ 

- -  x,6x 

i'.  support 


I 


Figure  16:  System  displacements  in  x.  Bodies  i  and  j  are  point  masses.  This  diagram 
is  used  to  determine  spring  forces  between  bodies  i  and  j  in  the  x  direction. 
The  flywheel  (j)  is  attached  to  the  end  of  a  rigid  link  length  d  extending 
from  the  support  structure  (i) 


are  assumed  to  be  small,  so  sin^  9  and  cos  9  ~  1.  x*  describes  the  displacement  of 
the  rotor  support  structure  in  the  x  direction.  The  rotor’s  position  in  x  is  dehned  as 
Xj.  The  X  direction  resultant  spring  force  on  mass  i  due  to  the  displacements  shown 
in  Figure  16  is  given  below  in  Equation  20. 


^Fx  T  9y^id  ^j) 


(20) 


Simlarly,  the  resultant  torques  on  mass  i  about  the  y  axis  are  described  below 
in  Equation  21. 


^Ty  ^y{^9y  i  ^yj')  kd(^9y  id  T  Xi  ^j) 


(21) 
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The  remaining  resultant  forces  and  torques  on  mass  i  can  be  determined  by 
analogy  or  by  the  construction  of  a  similar  diagram.  The  forces  and  torques  on  mass  i 
in  y  and  Ox  are  shown  below  in  Equation  22.  Since  this  model  ignores  rotor  motion  in 
2:  and  Oz^  those  equations  are  not  listed. 


TsFy  =  -kyijji  -  Ox, id  -  yj)  (22) 

^Tx  f^x{0x,i  dx,j^  kd(^9x^id  l/i  T  1/j) 

Likewise,  the  forces  acting  on  mass  j  (the  flywheel)  are  shown  below  in  Equa¬ 
tion  23. 


SFj,  =  -kx{xj  -Xi-  Oy^id)  (23) 

^y^Vj  Vi  Ox^id^ 


Equations  20-23  can  be  represented  in  matrix  form  as  shown  in  Equation  24 
where  i,  j  are  the  bodies  being  connected.  They  describe  the  forces  applied  to  the 
bodies  by  various  system  displacements. 


1 - 

1 

_  K,. 

- 1 

(24) 
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— kx  0  0  0  — kxd  0 

0  — kx  0  kxd  0  0 

0  0  0  0  0  0 

K,,,  = 

0  kxd  0  — kxd‘^  —  kx  0  0 

—kxd  0  0  0  —kxd?  —  kx  0 

0  0  0  0  0  0 

kx  0  0  0  0  0 

0  kx  0  0  0  0 

0  0  0  0  0  0 

K*,,  = 

0  —kxd  0  kt  0  0 

kxd  0  0  0  Kx  0 

0  0  0  0  0  0 


K  = 

— kx  0  0  0  0  0 

0  — kx  0  0  0  0 

0  0  0  0  0  0 

K,,,  = 

0  0  0  —  Kx  0  0 

0  0  0  0  —  Kx  0 

0  0  0  0  0  0 


For  this  model,  mass  i  is  always  the  support  structure,  so  i  =  1.  This  12  x  12 
stiffness  matrix  is  separated  into  four  parts  which  are  placed  into  the  global  system 
A  matrix  in  the  appropriate  locations. 
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The  damping  matrix,  C,  is  treated  the  same  way.  This  model  uses  damping 
proportional  to  stiffness  to  represent  the  losses  in  the  system.  This  is  represented  by 
c  =  eK,  where  ^  =  0.002,  which  is  reffected  by  the  values  in  Table  1. 

3.3.2  System  Equation  of  Motion  Assembly.  Recalling  Equations  2  and  17, 
a  state-space  system  EOM  is  written  as 

q  =  Aq  +  Bw 

Both  A  and  B  are  assembled  from  their  component  parts,  which  are  described 
in  Section  3.3.1.  This  can  also  be  represented  as  shown  in  Equations  25  and  26. 


A  =  (Ao  +  Ak  +  Ac) 


(25) 


B 


M 


-1 

system 


B 


system 


(26) 


where 


I  0 


M 


-1 

system 


0 


I  0 

0  m;' 
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0  I 
0  Gi 

Ag  = 

0  I 
0  G„ 


0  0  0  0 

SKi  1  0  Ki  2  0 

0  0  0  0 

Ak  K2,1  0  SK2j2  0 

0  0  0  0 

0  K„2  0 


0  0 

Kin  0 

0  0 

K2,n  0 

0  0 

SKn^n  0 


0  0  0  0 

0  SCi  1  0  Ci2 
0  0  0  0 

Ac  0  C2,l  0  SC2,2 

0  0  0  0 

0  Cn,l  0  Cn,2 
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Subscripts  on  stiffness  and  damping  terms  are  the  i,  j  subscripts  found  in  Equa¬ 
tion  24.  Placement  of  each  term  is  determined  from  Equation  18.  Ak  and  Ac  of 
Equation  25  are  more  general  than  the  model  requires.  Since  the  only  connections 
in  the  model  are  between  body  1  (the  support  structure)  and  other  bodies,  many  of 
the  terms  shown  in  these  general  equations  are  unnecessary,  leading  to  the  simplih- 
cation  of  Ak  as  Ak*  as  shown  below,  where  all  of  the  off-diagonal  terms  not  along 
the  first  row  or  column  are  zero.  Ac*  is  similar.  The  stiffness  term  for  the  support 
structure,  SKi^i,  also  contains  the  stiffness  between  the  fixed  satellite  bus  and  the 
support  structure. 


Ak*  = 


0 


0 


0  0  0 


SKii  0  K 


1,2 


K2,i  0  K 


0 


0  0  0 


2,2 


0 


0  0  0  0 

K„,1  0  0  0 


0  0 

Ki  „  0 

0  0 

0  0 

0  0 

K„,n  0 


The  input  matrix  B  is  simpler.  It  links  each  of  the  states  with  an  input  The 
column  and  row  of  a  B  term  determine  which  input  force  is  applied  to  which  state, 
respectively.  Because  the  EOM  is  represented  in  the  second  half  of  each  body’s  state 
vector,  that  is  where  the  inputs  are  applied. 
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B 


system 


0 

B„ 


(27) 


This  model  incorporates  rotating  centripetal  force  inputs  synchronized  to  the 
rotation  of  each  flywheel.  The  input  vector  describing  each  of  these  terms  is  shown  in 
Equation  28. 


u  = 


eica^  sin^it  eioj^  cos  6it 


sin  9nt  CnUJn  ^OS  1 


(28) 


As  Equation  28  shows,  the  input  vector  has  two  centripetal  force  terms  for  each 
wheel,  and  they  are  90°  out  of  phase  from  each  other.  The  magnitude  of  the  force  is 
Ciuf,  and  the  force  acts  in  the  direction  of  the  current  wheel  rotation,  di.  In  the  model 
this  is  divided  into  x  and  y  (or  9^  and  9y)  inputs,  which  vary  periodically  with  sin^j 
and  cos9i.  The  last  term,  1,  is  used  to  allow  for  a  constant  force  input  for  validation 
purposes.  This  vector  is  used  for  all  model  input.  Different  system  input  cases  are 
created  by  applying  forces  (imbalances  or  a  constant  force)  to  various  states  with 
changes  in  the  input  matrix  B. 


3.4  Integration 

The  model  created  in  Matlab®  is  integrated  numerically  using  ode45.  The 
differential  equation  created  by  the  state-space  model  is  kept  as  small  as  possible 
because  it  must  be  integrated  over  many  iterations.  Some  components  of  the  system 
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EOM,  however,  are  functions  of  time.  Both  the  gyroscopic  matrix  G  and  the  force 
vector  u  are  time  dependent  and  they  must  be  recomputed  for  each  iteration  of  the 
integration. 

In  addition,  careful  attention  must  be  paid  to  the  application  of  rotating  input 
forces.  A  periodic  rotating  input  force  is  represented  in  two  linear  dimensions  as  a  sine 
wave  in  one  dimension  and  a  cosine  wave  in  the  other  dimension.  If  these  forces  are 
directly  applied  to  an  unconstrained  mass  initially  at  rest,  they  can  cause  a  secular 
drift  in  one  direction.  This  phenomenon  occurs  in  the  model  used  for  this  thesis  (but 
not  the  actual  system).  One  solution  to  minimize  the  effect  of  the  secular  drift  is 
described  in  Appendix  A. 


3. 5  Additional  Components 

3.5.1  Appendage.  The  addition  of  a  flexible  appendage  is  accommodated  by 
adding  another  body  to  the  model.  An  appendage  is  used  to  study  the  beat  phe¬ 
nomenon  as  it  applies  to  this  model.  Since  the  studied  effects  of  this  appendage  are 
limited  to  9^  and  6y,  the  appendage  is  only  connected  to  the  model  with  torsional 
springs.  Also,  the  only  relevant  mass  properties  are  Ix  and  ly.  These  properties,  which 
were  chosen  to  give  the  appendage  a  natural  frequency  of  5  Hz,  are  shown  in  Table  3. 

Table  3:  Appendage  mass  properties 

mass  Ix,y  10  kg-m^ 

transverse  MOI  kt  9869  N-m/rad 


In  the  model,  the  appendage  is  added  as  a  new  body  with  mass  Mapp  = 


diag 


1  1  1  10  10  1 


and  stiffness  Kapp  dehned  according  to  Equation  24, 


with  all  terms  except  kt 
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3.5.2  Satellite/ support  spring.  When  required,  a  spring  is  connected  between 
a  firm  fixed  satellite  and  the  support  structure.  Since  input  forces  are  cyclical  and 
cause  primarily  rotations  in  the  support  structure,  it  is  the  rotational  DOFs  that  be¬ 


come  problematic  and  require  constraints.  This  additional  spring  is  attached  to  Ki^i, 


as  discussed  in  Section  3.3.2,  and  takes  the  form  Kgupport  =  — 1  x  10®diag 


0  0  0  1 


3. 6  Validation 

3.6.1  Validation  Inputs.  A  simple  IPACS  model  was  created  for  validation 
purposes.  The  flywheel  properties  for  this  model  are  shown  in  Table  4.  This  model 
has  a  long  rotor,  so  it  will  have  a  critical  coning  speed. 


Table  4:  Validation  model  flywheel  properties 


rotor  mass 

m 

10 

kg 

rotor  length 

1 

0.5 

m 

rotor  radius 

r 

0.12 

m 

support/rotor  distance 

d 

1 

m 

mag  bearing  stiffness 

k 

i^mag 

2500 

N/m 

polar  MOI 

Ip 

0.0781 

kg-m^ 

transverse  MOI 

It 

0.2474 

kg-m^ 

rotor  MOI  ratio 

P 

0.3158 

For  validation,  several  test  inputs  were  given  to  the  model  and  the  responses 
were  verified.  First,  a  few  constant-direction  forces  were  applied  to  ensure  that  signs 
were  correct  and  that  the  model  components  were  assembled  correctly.  Exponential 
rotational  growth  due  to  a  constant  applied  torque  is  shown  in  Figure  17.  Figure  17 
and  all  similar  figures  show  a  time  history  of  each  system  displacement  for  each  wheel. 
All  responses  are  shown  in  global  coordinates. 

The  next  set  of  tests  were  performed  to  ensure  the  model’s  consistency  in  ac¬ 
counting  for  wheel  spin  direction.  One  wheel  at  a  time  was  given  an  imbalance  input 
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(a)  Wheel  speed:  0 


(b)  Input  force:  as 
shown 


Figure  17:  System  demonstrating  exponential  displacement  growth  due  to  a  constant 
torque  input  applied  to  9x 
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■  vl  • 


to  test  the  response  of  that  wheel  and  the  system.  These  tests  are  shown  in  Figures  18 
and  19.  In  both  figures,  the  whole  system  is  rotating  together  slowly  as  a  result  of  the 
rotating  imbalance,  as  expected. 


X  io“'“  Support  wheel  1  ,( wheel  2 


(a)  Wheel  speed:  0/5  RPS  as  shown  (b)  Input  force:  as 

shown 

Figure  18:  System  demonstrating  forward  whirl 


(a)  Wheel  speed:  5/0  RPS  as  shown 


(b)  Input  force:  as 
shown 


Figure  19:  System  demonstrating  backward  whirl 


The  90°  phase  shift  between  the  Ox  and  6y  rotations  reveals  that  this  vibration  is 
a  coning  motion.  When  Ox  leads  Oy^  as  seen  in  Figure  18,  the  system  is  whirling  counter¬ 
clockwise  as  viewed  from  “above”  (looking  from  the  positive  z  direction  towards  the 
origin).  This  motion  is  a  forward  whirl  for  wheel  1,  since  that  rotor  is  spinning  in 
this  direction.  It  is  a  backward  whirl  for  wheel  2,  since  the  nominal  spin  direction  of 
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rotor  2  is  opposite  that  of  wheel  1.  When  6y  leads  O^i  the  system  is  rotating  in  the 
opposite  direction  and  the  whirl  orientations  are  reversed.  This  response  was  studied 
with  the  normal  system  model,  not  the  validation  model. 


Next  the  system’s  natural  frequencies  were  examined.  The  hrst  three  natural 
frequencies  of  a  rotor-bearing  system  are  given  in  Equations  11-13  and  15.  Recall  that 
fi  =  /2  are  bouncing  modes,  and  /s  and  /4  are  forward  and  backward  whirling  modes. 
The  system  responses  demonstrating  /1-/4  are  shown  in  Figures  20-24.  Figure  20 
demonstrates  the  bouncing  mode  fi  =  f2- 


(a)  Wheel  speed:  zero 


(b)  Initial  state: 
as  shown 


Figure  20:  System  demonstrating  bouncing  vibration  at  /i  ~  6.16  Hz 


Equation  11  gives  the  natural  frequency  of  a  spring- mass  system  as  fi  =  ^Jk/m/2^^. 
The  model’s  two  flywheels  were  displaced  symmetrically  to  make  the  system  equiv¬ 
alent  to  a  free  two-body  system.  The  equivalent  mass  of  a  freely  vibrating  two-body 
system  can  be  found  by  nieq  =  With  a  10  kg  support  structure  and 

two  10  kg  flywheels,  the  system  equivalent  mass  is  6.66  kg  and  /i  =  /2  ~  6.16  Hz. 
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Figure  21  shows  the  expected  forward  and  backward  whirl  speeds  of  the  valida¬ 
tion  model.  By  inspection,  when  wheel  speed  is  zero,  /s  =  /4  ~  5.66  Hz,  which  agrees 
with  Equation  12.  At  a  wheel  speed  of  10  RPS,  /s  ~  7.36  Hz  and  /4  ~  4.40  Hz. 


Figure  21:  Whirl  modes  of  the  validation  system.  Dotted  line  is  synchronous  with 
wheel  speed 

Figure  22  shows  that  the  forward  whirl  mode  behaves  as  expected  with  flywheel 
speed  a;  =  0.  Backward  whirl  at  this  speed  is  identical  except  for  direction,  which 
appears  as  an  opposite  phase  difference  in  6^  and  6y. 

The  whirl  modes  are  dependent  on  wheel  speed  and  they  diverge  with  increasing 
wheel  speed.  Figures  23  and  24  show  the  whirl  modes  when  wheel  speed  is  10  RPS. 
/3(10  Hz)  7.36  Hz  and  /4(10  Hz)  4.40  Hz.  In  both  cases  the  system  whirl  is  in 
the  same  direction  (counter-clockwise  as  viewed  from  “above”),  but  the  wheel  speed 
directions  were  changed  as  necessary  to  make  these  forward  or  backward  whirling 
motions. 
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(a)  Wheel  speed:  0  (b)  Initial  state:  as 

shown 


Figure  22:  Forward  whirl  with  zero  wheel  speed:  /a  =  /4  ~  5.66  Hz.  Backward  whirl 
is  identical  except  for  direction 
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(a)  Wheel  speed:  10  RPS  as  shown 


(b)  Initial  state:  as 
shown 


Figure  23:  Forward  whirl  mode  at  /s  ~  7.36  Hz  when  wheel  speed  is  10  RPS 
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(a)  Wheel  speed:  10  RPS  as  shown 


(b)  Initial  state:  as 
shown 


Figure  24:  Backward  whirl  mode  at  /4  ~  4.40  Hz  when  wheel  speed  is  10  RPS 


Forward  whirl  is  speed  dependent  and,  for  a  long  rotor,  there  is  a  critical  con¬ 
ing  speed  where  the  whirling  frequency  is  the  same  as  the  wheel  frequency.  Figure  25 
demonstrates  this  phenomenon.  Wheel  speed  is  varied  from  4-12  RPS,  and  vibration 
is  much  worse  at  the  critical  coning  frequency  fcon  ~  6.8  Hz.  The  difference  between 
the  support  structure’s  vibration  in  6x  and  6y  is  a  modeling  artifact. 


3. 7  Summary 

The  model  developed  for  this  thesis  was  tested  and  performed  correctly  given  a 
variety  of  input  scenarios,  including  initial  conditions,  constant  forces,  and  periodic 
forces.  Critical  coning  frequency  was  demonstrated  successfully,  as  well  as  forward 
and  backward  whirls. 
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Figure  25:  Model  demonstrating  the  critical  coning  frequency.  As  the  wheel  speed 
passes  through  fcon  ~  6.8  Hz,  the  vibration  amplitude  gets  much  larger. 
Wheel  speed  varies  linearly  from  4  RPS  at  t=4  to  12  RPS  at  t=12 
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IV.  Results  and  Analysis 

Several  tests  were  conducted  using  the  model  described  in  this  thesis,  and  the  results 
are  presented  here.  First,  a  full  frequency  sweep  was  performed  to  see  if  any  imme¬ 
diate  problems  arose.  Other  tests  examined  the  existence  of  beat  phenomenon  and 
extra-synchronous  whirl  excitation.  For  this  section,  flywheel-specihc  rotation  speed 
is  represented  as  Un- 


4-1  Full  Envelope  Sweep 

First  a  variety  of  full  envelope  frequency  sweep  tests  were  carried  out.  For  these 
tests,  wheel  speed  of  one  or  both  rotors  varied  from  333-1000  RPS.  A  moment- 
inducing  imbalance  in  each  rotor  caused  a  vibration  input  synchronized  with  each 
wheel  speed. 

These  tests  revealed  nothing  of  interest  except  that  the  system’s  gyroscopic 
stiffness  diminishes  as  the  wheel  speeds  approach  equality  in  opposite  directions.  This 
reduction  in  gyroscopic  stiffness  is  illustrated  nicely  in  Figures  26-28. 


X  .,0-^°  Support 


X  1 0”''°  wheel  1  x  1 0"^°  wheel  2 


X  10 


-10 


X  10 


-10 


Figure  26:  System  rotations  as  a  result  of  moment-inducing  imbalance  vibrations. 
Wheel  speeds:  uii  =  333,  U2  =  333-1000  RPS 
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Figure  27:  System  rotations  as  a  result  of  moment-inducing  imbalance  vibrations. 
Wheel  speeds:  Wheel  speeds:  cui  =  1000,  002  =  333-1000  RPS 
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Figure  28:  System  rotations  as  a  result  of  moment-inducing  imbalance  vibrations. 
Wheel  speeds:  cui  =  1000-333,  002  =  333-1000  RPS 
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The  system  response  clearly  changes  frequency  in  each  case.  The  time-varying 
long-period  wobble  in  the  system  is  caused  by  the  sudden  application  of  the  input 
forces,  which  causes  a  slight  misalignment  between  the  angular  momentum  vector 
and  the  flywheel  rotation  axes.  The  system  angular  momentum  vector  gets  larger  as 
the  flywheel  speeds  get  further  apart,  giving  the  system  as  a  whole  more  gyroscopic 
stiffness.  The  additional  stiffness  manifests  in  an  increased  oscillatory  frequency  and 
a  smaller  displacement  from  the  equilibrium  state  as  expected. 

In  all  cases,  the  vibration-inducing  imbalance  forces  are  less  signihcant  than  the 
system’s  slow  oscillation.  The  magnitude  of  the  imbalance-induced  vibration  is  evi¬ 
denced  by  the  thickness  of  the  lines  seen  in  Figures  26-28.  A  closer  look  at  the  rotation 
of  the  support  structure  as  shown  in  Figure  29  reveals  that  the  wheel-synchronous 
vibration  is  small  relative  to  the  long-period  oscillation.  This  hgure  shows  the  same 
data  seen  in  Figure  26. 

Adding  a  flexible  appendage  (representing  a  satellite  structure  such  as  an  an¬ 
tenna  or  solar  array)  to  the  system  does  not  change  the  response  very  much.  Figure  30 
shows  the  response  of  the  system’s  flexible  appendage,  which  has  a  natural  frequency 
of  5  Hz.  There  is  some  vibration  evident  at  the  natural  frequency  of  the  appendage, 
but  the  response  is  dominated  by  the  slow  oscillation  of  the  support  structure. 

The  slow,  changing  oscillation  seen  in  Figures  26-30  accurately  reflects  the  re¬ 
sponse  of  a  system  with  changing  gyroscopic  stiffness  given  an  initial  whirling  type 
motion.  The  initial  whirling  motion  seen  in  these  hgures,  however,  is  an  artifact  of 
modeling,  caused  by  the  instant  application  of  a  disturbing  force.  Such  a  motion  is 
possible  and  indeed  uncontrollable  in  a  two  wheel  IPACS,  but  it  can  only  be  caused 
by  an  external  disturbance  force.  In  an  IPACS  with  full  3-axis  control  authority,  the 
control  in  other  dimensions  would  be  used  to  eliminate  the  whirling  motion. 
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X  10“" 


Figure  29:  A  closer  look  at  the  support  structure  vibration  from  Figure  26.  This  view 
shows  the  wheel-synchronous  vibration 


X  10 


time  (s) 

Figure  30:  A  closer  look  at  the  appendage  vibration  from  Figure  26.  This  view  reveals 
a  small  amount  of  vibration  at  the  appendage’s  natural  frequency,  but 
behavior  is  dominated  by  long-period  oscillation 
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4-2  Beat  Frequency  Analysis 


After  studying  the  effect  of  near-frequency  input  vibrations,  it  is  evident  that 
this  model  does  not  reveal  any  beat  frequency  problems.  This  is  unsurprising  given 
the  linear  nature  of  the  model. 

The  flexible  appendage  did  experience  a  small  resonant  vibration  when  the 
difference  of  two  input  frequencies  matched  the  natural  frequency  of  the  appendage, 
but  it  experienced  the  same  vibration  with  other  unrelated  input  frequencies. 

First,  the  two-wheel  model  was  run  with  slightly  different  wheel  speeds:  333  and 
338  RPS.  The  5  RPS  difference  in  the  wheel  speeds  creates  a  beat  frequency  which 
matches  the  natural  frequency  of  the  flexible  appendage,  which  was  not  yet  added. 
As  Figure  31  shows,  the  5  Hz  beat  frequency  is  clearly  evident. 
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Figure  31:  Beat  frequency  is  clearly  visible  in  the  support  structure  when  oji  = 
333,  UJ2  =  338  RPS 

Next  the  model  was  run  at  the  same  wheel  speeds  with  the  flexible  appendage 
added.  The  appendage  represents  a  flexible  satellite  body  such  as  a  solar  panel  or 
antenna.  The  rotational  state  history  from  this  test  is  shown  in  Figure  32.  The  beat 
frequency  is  still  clearly  visible,  and  the  appendage  is  experiencing  vibration  near  its 
natural  frequency. 
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Figure  32:  Beat  frequency  is  still  visible  in  the  support  structure  when  oji  =  333,  0J2  = 
338  RPS,  and  the  appendage  experiences  excitation  at  its  resonant  fre¬ 
quency  of  5  Hz 


Finally,  the  same  model  with  flexible  appendages  was  run  with  a  new  and  wider 
wheel  speed  difference:  uji  =  333,  002  =  400  RPS.  This  is  shown  in  Figure  33.  The 
beat  frequency  is  no  longer  present,  but  the  appendage  still  exhibits  vibration  near  its 
natural  frequency.  The  vibration  magnitude  of  the  appendage  in  this  hgure  is  largely 
the  same  as  seen  in  Figure  32. 


time  (s) 


Figure  33:  Beat  frequency  is  gone,  but  appendage  still  vibrates  near  its  resonant  fre¬ 
quency  of  5  Hz  when  ui  =  333,  U2  =  400  RPS 


From  studying  the  results  of  the  tests  performed  with  both  frequency  offsets, 
it  appears  that  there  is  no  more  low  frequency  excitation  at  the  difference  frequency 
of  two  inputs  than  there  is  at  any  other  combination  of  inputs.  The  beat  frequency 
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exists,  but  its  effects  are  very  limited.  Recalling  Equation  16,  the  combined  output  of 
two  signals  that  are  close  in  frequency  can  be  written  as 


y  =  2  sin 


t  cos 


t 


The  maximum  magnitude  of  the  resultant  vibrations  shown  in  Figure  31  should 
be  twice  that  of  the  vibrations  evident  when  only  one  source  of  vibration  is  present. 
To  check  this,  a  test  was  run  with  only  one  wheel  spinning.  Figure  34  shows  that 
the  vibration  magnitude  is  indeed  approximately  half  of  the  peak  seen  when  a  beat 
frequency  is  present. 


time  (s) 


Figure  34:  Vibration  here,  caused  by  one  wheel  rotating  at  333  RPS,  is  approximately 
half  the  peak  vibration  seen  in  Figure  31 


A  power  spectral  density  plot  of  the  support  structure’s  vibrations  during  beat 
frequency  and  non-beat  frequency  conditions  reveals  that  there  is  no  power  at  the 
low  beat  frequency.  During  the  beat  case,  wheel  speeds  are:  ui  =  333,  UJ2  =  338  RPS. 
For  the  non-beat  case,  the  speed  of  wheel  2  is  raised  to  UJ2  =  400  RPS.  The  peak- 
normalized  PSD  of  both  responses  is  shown  in  Figure  35. 

The  PSD  reveals  that  there  is  no  energy  at  the  low  beat  frequency,  and  the 
rotation  time  histories  reveal  that  there  is  no  extraneous  excitation  as  a  result  of 
the  beat  frequency.  Aside  from  periodically  doubling  the  magnitude  of  the  resultant 
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Figure  35:  Peak-normalized  PSD  of  support  structure  9^  during  beat  (cui  =  333,  U2  = 
338  RPS)  and  non-beat  (cui  =  333,  U2  =  400  RPS)  conditions.  Low  fre¬ 
quency  response  is  nearly  identical.  Note  peaks  at  the  input  excitation 
frequencies 
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high-frequency  vibration,  there  is  no  negative  effect  from  the  beat  phenomenon  in  a 
linear  system. 

4.3  Extra- Synchronous  Whirl  Excitation 

Long  rotors  have  a  critical  coning  frequency,  fcon  where  the  wheel  speed  is 
synchronous  with  the  forward  whirl  speed.  This  was  illustrated  with  Figure  25  in  Sec¬ 
tion  3.6.1.  All  rotors,  however,  have  forward  and  backward  whirl  modes.  With  multi¬ 
ple  flywheels  in  one  system,  the  flywheels  provide  each  other  with  extra-synchronous 
vibration  sources,  and  it  may  be  possible  for  these  extra-synchronous  frequencies 
to  excite  unexpected  resonances.  Both  sub-  and  super-synchronous  excitations  were 
studied  for  this  phenomenon. 

In  this  section,  uj  refers  to  flywheel-specific  wheel  speed.  Forward  and  backward 
whirl  speed  are  respectively  referred  to  as  /a  and  /4,  which  were  derived  and  explained 
in  Section  2.3.5.  Spectrograms  are  used  in  this  section  to  illustrate  the  time-varying 
nature  of  the  response  frequencies.  All  spectrograms  are  derived  from  9^  (the  rotation 
about  x)  of  the  wheel  in  qnestion.  Units  are  very  small  and  are  included  only  for  the 
purpose  of  comparison  between  different  cases.  The  satellite/support  spring  was  used 
for  all  of  the  tests  in  this  section. 

4.3.1  Sub- Synchronous  Whirl  Excitation.  A  long-rotor  model  was  used  to 
look  for  sub- synchronous  whirl  excitation.  Sub-synchronous  resonance  is  vibration 
of  a  long  rotor  at  speeds  lower  than  the  wheel  speed.  Snb-synchronous  resonance 
is  at  a  higher  frequency  than  the  critical  coning  frequency,  /con,  where  the  forward 
whirl  frequency  crosses  the  wheel  frequency.  The  region  of  interest  using  the  flywheel 
parameters  chosen  for  this  study  is  seen  in  Figure  36,  which  shows  both  the  forward 
and  backward  whirl  of  the  long  rotor  used  in  this  study. 


Figure  36:  Forward  and  backward  whirl  modes  of  the  sub-synchronous  study  model 

In  the  sub-synchronous  test,  the  wheels  will  spin  at  two  different  speeds.  The 
vibration  input  of  the  slower  wheel  should  excite  a  forward  whirling  mode  in  the  faster 
wheel,  even  if  the  faster  wheel  is  perfectly  balanced. 

The  flywheel  parameters  used  for  this  study  are  shown  in  Table  5.  Other  param¬ 
eters  remain  the  same  as  those  shown  in  Table  1  from  Section  3.2.3.  This  model  was 
studied  both  with  and  without  bearing  damping  included.  These  values  were  chosen 
not  for  realism,  but  to  give  favorable  conditions  for  excited  whirl  at  a  frequency  other 
than  the  wheel  frequency. 

Table  5:  Flywheel  parameters  for  sub-synchronous  study 


mass 

m 

10 

kg 

polar  MOI 

Ip 

0.1125 

kg-m' 

transverse  MOI 

It 

0.2015 

kg-m’ 

rotor  MOI  ratio 

P 

0.5583 
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In  the  sub- synchronous  test,  wheel  2  was  perfectly  balanced  and  operated  at 
700  RPS.  This  gives  it  a  forward  whirl  frequency  of  /s  ~  400  Hz  and  a  backward 
whirl  frequency  of  ~  —11  Hz.  Wheel  1  was  given  a  moment-inducing  imbalance 
and  operated  at  speeds  intended  to  excite  /s  of  the  faster  balanced  wheel,  or  about 
400  RPS.  Wheel  I’s  own  speed,  ui  caused  separate  whirling  frequencies  in  wheel  1: 
/g  Ri  242  Hz  and  /4  ~  —18  Hz.  Natural  frequencies  can  be  found  by  inspection  from 
Figure  36  and  analytically  with  Equations  12  and  13  from  Section  3.2.3.  In  all  cases 
the  speed  of  wheel  1  (the  unbalanced  excitation  wheel)  varies  from  350  RPS  at  t  =  0 
to  450  RPS  at  f  =  2. 


4-3. 1.1  Undamped  sub- synchronous  response.  At  hrst,  damping  was 
neglected  to  see  if  there  would  be  any  response.  Wheel  I’s  response  is  shown  in 
Figure  37.  Note  the  strength  of  the  self-induced,  sub-synchronous,  forward  whirling 
mode.  This  wheel  is  well  past  its  critical  coning  speed  for  forward  whirl,  but  it  is 
still  responding  strongly  at  the  whirling  frequency.  In  fact,  at  lower  wheel  speeds  the 
resonant  forward  whirling  response  is  stronger  than  the  input  disturbance. 

-23 

x10 


Figure  37:  Spectrogram  of  wheel  I’s  response  to  its  own  undamped  imbalance,  oji 
varies  from  350-450  RPS.  Note  the  strong  forward  whirling  resonance 
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Wheel  2’s  response  over  the  same  time  period  was  examined  at  multiple  speeds. 
Figure  38  shows  the  wheel  1  response  at  a  baseline  of  wi  =  0  and  at  ca  =  =f700  RPS. 
Wheel  2  experienced  a  signihcant  forward  whirling  motion  when  it  was  rotating  in 
the  same  direction  as  wheel  1.  Support  response  in  all  cases  was  similar. 

Two  spurious  vibrations  can  be  found  in  Figure  38.  Most  notable  is  the  constant 
frequency  vibration  at  approximately  50  Hz  which  appears  in  38(a)  and  38(c).  The 
50  Hz  vibration  appears  to  be  worse  when  the  excitation  is  at  400  Hz,  but  it  is  actually 
worse  at  approximately  1  sec,  and  occurs  repeatedly  if  the  integration  time  is  longer. 
There  is  also  a  vibration  at  approximately  70  Hz  visible  in  38(a).  These  vibrations  do 
not  match  either  of  the  rotor’s  resonant  frequencies,  and  they  are  system  level  effects 
caused  by  the  bus/support  springs  and  the  net  gyroscopic  stiffness. 

Two  notable  vibrations  do  appear,  however.  The  first  of  these  vibrations  is  at 
the  backward  whirl  frequency,  /4,  of  wheel  1  (the  excitation  wheel).  This  vibration  at 
approximately  20  Hz  is  visible  in  all  three  scenarios  of  Figure  38,  but  is  strongest  in 
the  response  shown  in  Figure  38(c),  which  is  the  case  where  the  wheels  are  spinning  in 
opposite  directions.  This  excitation  indicates  that  wheel  1  is  exciting  its  own  backward 
whirl  mode  through  interaction  with  wheel  2,  even  though  wheel  2  is  perfectly  balanced. 

The  other  vibration  apparent  in  Figure  38  is  the  vibration  that  appears  near 
t  =  1  second  at  approximately  400  Hz  in  38(b).  This  frequency  is  the  forward  whirl 
speed  of  wheel  2,  and  the  intended  excitation  frequency  of  this  test  model.  This  vibra¬ 
tion  appears  and  is  most  severe  about  halfway  through  the  test,  when  the  excitation 
frequency  is  400  Hz.  Forward  whirl  in  wheel  2  is  expected  to  be  worst  at  that  time 
because  it  is  being  excited  at  /a.  The  persistence  of  the  whirling  motion  after  the 
excitation  frequency  has  moved  past  the  resonant  frequency  shows  that  whirl  tends 
to  persist  after  excitation.  The  same  time  history  appeared  when  wheel  I’s  speed  was 
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Figure  38:  Spectrogram  of  wheel  2’s  response  to  sub-synchronous  vibration  input. 

uji  varies  from  350-450  RPS.  Note  the  different  frequency  scale  in  subhg- 
ure  (b).  Vibration  at  50  Hz  fades  in  and  out  repeatedly 
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varied  backwards,  from  450-350  RPS.  In  other  words,  the  appearance  and  decrescendo 
of  a  whirling  motion  was  consistent  in  time,  not  in  input  frequency. 

Note  that  the  model  conhguration  represented  in  38(b)  is  not  operationally 
realistic  for  a  two  wheel  IPACS,  but  it  demonstrates  sub-synchronous  excitation.  The 
problem  is  that  both  wheels  are  spinning  in  the  same  direction,  which  does  not  allow 
for  high  power  storage  and  transfer  with  a  net  zero  angular  momentum  change.  In 
reality,  the  wheels  will  be  spinning  in  opposite  directions  as  they  are  for  the  case 
shown  in  Figure  38(c).  This  situation  might  occur,  however,  in  a  bank  of  multiple 
rotors,  such  as  one  being  used  primarily  for  energy  storage. 

4 .3. 1.2  Damped  sub-synchronous  response.  Since  sub-synchronous  ex¬ 
citation  was  found  to  affect  a  two  flywheel  system  with  no  damping,  the  next  step 
was  to  see  if  it  remained  when  damping  was  included  in  the  model.  The  same  test 
conditions  used  in  the  undamped  case  were  used  to  study  a  more  realistic  damped 
flywheel  conhguration. 

Figure  39  shows  the  response  of  wheel  1,  the  excitation  wheel.  Note  that  the 
resonant  whirling  response  has  been  attenuated  below  the  threshold  of  visibility  in 
this  spectrogram  and  all  that  remains  is  the  synchronous  imbalance-induced  vibration. 
Since  the  wheel  speed  is  so  far  from  the  critical  whirling  frequency  fcon,  even  the  very 
small  amount  of  damping  present  in  this  model  is  enough  to  eliminate  this  resonance. 

Figure  40  shows  the  damped  response  of  wheel  2  to  the  same  sub-synchronous 
imbalance  input  that  induced  forward  and  backward  whirl  in  the  undamped  case.  All 
resonances  have  fallen  below  the  detectable  threshold,  leaving  only  the  direct  input 
frequency  created  by  the  unbalanced  rotation  of  wheel  1. 

From  the  signihcant  and  benehcial  effect  of  even  modest  damping  in  this  multiple 
flywheel  IPACS  model,  it  appears  that  real  flywheel  systems  are  safe  from  the  harmful 
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Figure  39:  Spectrogram  of  wheel  I’s  response  to  its  own  damped  imbalance,  uji  varies 
from  350-450  RPS.  The  forward  whirling  resonance  is  gone 

effects  of  sub-synchronous  flywheel  vibration  in  an  axially  opposed  conhguration. 
Real  flywheels  necessarily  include  damping  since  it  is  a  real  and  unavoidable  physical 
phenomenon. 


4-3.2  Super- synchronous  Whirl  Excitation.  A  short-rotor  model  with  the 
nominal  parameters  given  in  Table  1  from  Section  3.2.3  was  used  to  look  for  super- 
synchronous  whirl  excitation.  Super-synchronous  resonance  is  vibration  of  a  short 
rotor  at  speeds  higher  than  the  wheel  speed.  Ordinarily  this  is  not  a  problem  for  short 
rotors,  since  their  forward  whirling  frequencies  always  remain  higher  than  the  wheel 
speed.  An  external  vibration  input  at  the  resonant  speed  could  induce  vibration  in 
this  mode,  however.  Figure  41  shows  the  forward  and  backward  whirling  frequencies, 
/a  and  f^,  for  this  configuration. 

This  test  used  the  same  setup  as  the  test  for  sub-synchronous  whirl,  with  dif¬ 
ferent  rotors  and  a  different  excitation  speed.  In  this  test  wheel  1  was  spun  from 
950-1050  RPS  in  an  attempt  to  excite  vibration  modes  in  wheel  2,  which  was  spin¬ 
ning  at  =f700  RPS.  A  non-spinning  wheel  was  again  used  as  a  baseline.  Whirl  modes 
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Figure  40:  Wheel  2  damped  sub-synchronous  spectrograms,  oji  varies  in  time  from 
350-450  RPS.  Frequencies  other  than  the  input  have  been  attenuated 
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Figure  41:  Forward  and  backward  whirl  modes  of  the  super-synchronous  study  model. 

One  dotted  line  is  zero  and  the  other  is  synchronous  with  the  wheel  speed. 

for  wheel  2  are  /a  ~  890  Hz  and  fi  ~  —12  Hz.  For  wheel  1,  /4  ~  —7.9  Hz.  /s  for 
wheel  1  is  well  above  any  vibrations  present  in  the  system,  since  the  rotors  are  short. 

Damping  was  again  at  first  neglected  to  see  if  there  would  be  a  response.  The 
super-synchronous  excitation  response  was  very  similar  to  the  sub-synchronous  excita¬ 
tion  response.  Spectrograms  in  this  section  show  a  wheel  1  speed  range  from  850  RPS 
at  f  =  0  to  950  RPS  at  f  =  2. 

Spectrograms  of  the  wheel  1  and  support  responses  did  not  reveal  any  reso¬ 
nances.  The  only  place  in  the  model  with  any  visible  extra-synchronous  frequency 
response  was  wheel  2.  The  spectrograms  for  wheel  2  are  shown  in  Figure  42.  Simi¬ 
lar  to  the  sub-synchronous  case,  there  is  a  slight  forward  whirl  response  when  spin 
directions  are  identical,  and  there  is  a  slight  backward  whirl  excitation  regardless  of 
relative  spin  directions. 
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Figure  42:  Spectrograms  of  wheel  2’s  undamped  response  to  super-synchronous  vi¬ 
bration  input.  Ui  varies  in  time  from  850-950  RPS 
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Figure  42  also  reveals  a  slight  response  around  the  backward  whirl  modes.  There 
is  not  enough  spectral  resolution  in  the  hgures  to  determine  which  wheel’s  response 
is  causing  the  vibration.  They  are  very  close  in  frequency,  however  (approximately  8 
and  12  Hz),  and  perhaps  both  of  them  are  responding. 

As  seen  in  Figure  43,  resonant  vibrations  due  to  a  super-synchronous  vibration 
excitation  disappear  when  damping  is  included  in  the  model,  similar  to  the  disap¬ 
pearance  of  resonant  vibration  in  the  sub-synchronous  test. 

4-4  Summary 

The  model  developed  in  this  thesis  was  used  to  investigate  potential  causes  of 
troublesome  low  frequency  vibration  in  an  IPACS.  Beat  frequency  was  found  to  have 
absolutely  no  impact  on  low  frequency  vibrations,  and  minimal  impact  at  any  fre¬ 
quency,  with  the  only  impact  being  the  periodic  doubling  of  the  input  forces.  Extra- 
synchronous  resonance  excitation,  both  sub-  and  super-synchronous,  was  found  to 
have  only  a  small  impact  on  satellite  vibrations.  Backward  whirl  was  excited  regard¬ 
less  of  spin  directions,  but  extra-synchronous  forward  whirl  was  only  excited  when 
spin  directions  were  identical.  Extra-synchronous  vibration  is  very  small  relative  to 
the  inducing  imbalance,  but  it  is  at  a  lower  frequency.  The  identical  spin  direction 
conhguration  is  of  no  concern  for  a  simple  IPACS,  but  could  be  a  factor  in  a  bank 
of  flywheels  intended  primarily  for  energy  storage.  With  a  small  amount  of  damping, 
however,  all  extra-synchronous  vibration  effects  disappear. 
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Figure  43:  Spectrograms  of  wheel  2’s  damped  response  to  super-synchronous  vibration 
input,  uji  varies  in  time  from  850-950  RPS 
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V.  Conclusion 


5. 1  Summary 

In  this  thesis  a  flexible  linear  dynamics  model  was  developed  to  study  vibra¬ 
tions  in  a  system  of  multiple  connected  gyroscopic  flywheels.  The  model  was  validated 
through  several  types  of  test  input,  and  then  used  to  analyze  various  sources  of  vibra¬ 
tion  interaction  between  two  flywheels  a  simulated  IPACS.  System  parameters  were 
intended  to  be  representative  of  a  near-future  technology  demonstration  mission  to 
validate  the  use  of  flywheels  for  energy  storage  and  attitude  control  in  space. 

The  source  of  vibration  studied  in  this  model  was  an  assumed  eccentricity  caused 
by  manufacturing  defect.  Furthermore,  it  was  assumed  that  this  defect  was  attenuated 
by  active  magnetic  bearing  control  methods  described  in  previous  literature.  The 
residual  vibrations  after  attenuation  were  used  as  model  inputs.  The  vibration  inputs 
were  applied  in  the  radial  and  transverse  directions  to  an  axially-aligned  two-flywheel 
IPACS. 

The  model  was  then  used  to  study  interactions  between  multiple  sources  of 
vibration  in  an  IPACS.  Specifically,  the  beat  frequency  and  extra-synchronous  vibra¬ 
tion  excitation  phenomena  were  studied,  including  both  sub-  and  super-synchronous 
vibration  excitation  scenarios. 

5. 2  Findings 

The  beat  frequency  problem  was  shown  to  have  almost  no  impact  in  a  linear 
system.  The  extra-synchronous  vibration  effects  were  both  found  to  have  a  small  effect 
on  a  double  flywheel  system.  When  mild  damping  was  applied,  effects  were  below  the 
threshold  of  observability,  but  with  no  damping  there  was  a  small  but  consistent 
system  response  to  backward  whirl  modes.  In  addition,  if  there  are  multiple  flywheels 
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with  the  same  spin  direction,  extra-synchronous  forward  whirl  in  one  wheel  can  be 
excited  by  another  unbalanced  wheel. 

Of  the  vibration  sources  considered  in  this  thesis,  extra-synchronous  resonance 
excitation,  both  sub-  and  super-synchronous,  proved  to  be  of  some  potential  concern 
in  the  future  development  of  IPACS.  As  bearing  and  flywheel  system  technology  con¬ 
tinues  to  improve,  damping  will  continue  to  fall  as  flywheel  designers  seek  to  remove 
sources  of  loss  from  their  systems.  At  some  point  the  damping  may  be  small  enough 
to  allow  the  extra-synchronous  resonances  to  affect  the  system  negatively.  Vibrations 
caused  by  input  interactions  are  much  smaller  than  the  vibrations  themselves,  but 
they  are  also  at  a  lower  frequency. 

Park  proposed  a  notch  hlter  centered  at  the  wheel  speed  to  mitigate  the  effects  of 
wheel  imbalance  to  an  IPACS  or  satellite.  If  bearing  development  produces  magnetic 
bearings  with  low  enough  damping  to  allow  the  extra-synchronous  resonant  vibration 
to  become  a  problem,  the  control  system  designer  could  also  develop  methods  for 
mitigating  extra-synchronous  vibration  modes. 

5.3  Contributions 

This  thesis  proved  the  existence  of  inter-flywheel  vibration  interactions  in  a 
multiple  flywheel  system.  Extra-synchronous  resonance  excitation  between  multiple 
rotors  was  found  to  exist  in  an  ideal  undamped  conhguration,  but  even  a  very  small 
realistic  amount  of  damping  was  enough  to  mitigate  the  effect  to  the  point  that  it 
was  less  of  a  concern  than  individual  rotor  vibration  inputs. 

5.4  Recommendations  for  Future  Work 

This  thesis  showed  that  damping  can  mitigate  and  possibly  eliminate  unde¬ 
sirable  extra-synchronous  flywheel  vibration.  Future  research  should  be  done  to  de- 
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termine  how  low  the  bearing  damping  can  be  before  the  sub-synchronous  resonance 
negatively  affects  system  operation. 

In  addition,  the  model  developed  in  this  thesis  was  used  to  investigate  some  of 
the  effects  of  multiple  flywheel  vibrations  in  an  IPACS,  but  this  thesis  was  limited 
in  scope  to  analyzing  two  co-axial  rotors  with  simple  imbalance-induced  vibrations 
in  and  about  2  axes.  In  an  advanced  flywheel  there  is  also  signihcant  potential  for 
rotational  vibration  about  the  spin  axis  due  to  torque  ripple.  In  the  co-axial  case  the 
torsional  vibrations  will  remain  independent  from  other  system  disturbances,  but  in 
a  flywheel  system  with  multiple  non-co-axial  rotors,  torsional  vibrations  will  induce 
other  vibrations  that  may  be  deleterious  to  successful  IPACS  implementation  on  a 
satellite. 

The  analysis  of  those  extra  factors  should  be  undertaken  as  a  further  step  along 
the  path  to  developing  flywheel  energy  storage  systems.  With  torsional  stiffness  turned 
on  and  the  appropriate  input  forces  applied,  the  model  developed  in  this  thesis  is 
capable  of  performing  this  evaluation. 

Furthermore,  as  flywheel  performance  continues  to  improve  and  bearing  de¬ 
velopment  produces  magnetic  bearings  with  lower  losses,  the  vibration  interaction 
effects  uncovered  in  this  thesis  may  become  signihcant.  If  that  occurs,  hywheel  bear¬ 
ing  controllers  must  be  able  to  account  for  and  mitigate  the  effects  of  backward  whirl 
excitation  between  hywheels. 

5. 5  Conclusion 

Flywheels  may  someday  be  a  key  component  for  providing  the  satellite  designer 
with  hexibility  and  performance  benehts.  Their  successful  on-orbit  deployment  has 
to  this  point  been  delayed  longer  than  any  of  the  original  proponents  of  the  idea 
could  have  envisioned,  and  many  problems  in  this  area  remain  unstudied.  Someday, 
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however,  enough  problems  will  be  solved  to  make  flywheels  in  orbit  a  viable  reality. 
The  engineering  community  must  continue  to  work  diligently  in  this  area  until  that 
day  arrives. 
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Appendix  A.  One  Integration  Problem  and  a  Solution 

Position  is  the  double  integration  of  acceleration,  which  is  proportional  to  the  input 
force  as  shown  in  Equation  29. 

x{t)  =  J  v(t)  ^  jj  ait)  =  ^  (29) 

Given  sinusoidal  inputs  acting  on  a  1  kg  mass: 


fit)  = 

sint 

cost 

a{t)  = 

sint 

cost 

v{t)  = 

—  cos  t  +  Cl 

sin  t  +  Cl 

xit)  = 

—  sint  +  cit  +  C2 

—  cost  +  Cit  +  C2 

If  the  mass  is  initially  at  rest  (xq  =  Uq  =  Oq  =  0),  the  constants  are  found  to  be: 


fit)  = 

sint 

cost 

Cl  = 

1 

0 

C2  = 

0 

1 

When  a  sine  wave  force  is  applied  to  an  unconstrained  mass  initially  at  rest, 
Cl  will  cause  a  positive  secular  drift.  This  drift  does  not  occur  if  the  input  force  is 
a  cosine  wave.  The  difference  between  sine  and  cosine  input  forces  is  illustrated  in 
Figure  44. 

The  drift  illustrated  in  Figure  44  is  a  phenomenon  observed  only  in  modeling; 
it  does  not  occur  in  real  systems.  Most  models  do  not  reveal  the  drift  because  they 
model  systems  constrained  by  some  spring  force.  For  unconstrained  systems,  a  model 
can  reveal  a  drift,  but  it  does  not  exist  in  reality  because  the  periodic  rotating  forces 
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Figure  44:  Acceleration,  velocity,  and  position  due  to  sinusoidal  inputs  applied  to  a 
free  body  initially  at  rest 

such  as  wheel  imbalances  always  grow  from  an  initial  magnitude  of  zero.  It  is  only  a 
factor  when  modeling  a  system  that  is  rotating  at  t  =  0. 

Because  the  rotating  force  must  be  represented  as  a  sine  wave  in  one  dimension, 
the  direct  application  of  this  force  will  cause  a  secular  drift.  One  solution  is  to  delay 
the  application  of  this  force  until  the  input  force  has  rotated  90°.  This  delay  will  cause 
the  sine  wave  to  look  like  a  cosine  wave  (it  will  begin  from  1  rather  than  from  0), 
and  it  will  be  shifted  correctly  in  phase.  With  multiple  input  forces,  each  one  must 
be  “turned  on”  at  the  appropriate  moment.  Given  the  numerical  nature  of  the  model 
in  this  thesis,  the  limits  of  integration  must  match  these  times  as  well  to  ensure  that 
the  sine  input  is  applied  at  exactly  the  right  moment.  “Turning  on”  the  forces  at  the 
correct  times  will  eliminate  this  source  of  drift. 
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Appendix  B.  Matlab®  Code 

Listing  B.l;  ../models/model.m 

function  [output]  =  model  (config) 
global  ss  model  input  output  c  ii  model 
7,  Feb  2011 
7«  Jordan  Firth 

7«  state  =  [x  y  z  a  b  g  x’  z’  a’  b’  g’]^ 

7«7«  constants 

7«  compute/import  data  from  case_maker 
[model , input ]  =  case_r eader ( conf ig ) ; 

7«7«  R  Matrices 

7«  r3-  each  page  of  array  is  rot  mat  for  one  wheel 

7«  rl2-  blkdiag  of  4  r*3’s-  rotate  all  4  states  of  the  wheel  at  once 

7«  R“  one  giant  rotation  matrix  for  the  whole  massive  A  matrix-may  not  need 

7«  Here  are  the  rz  ’  s 

r3  =  zeros  (3 , 3  ,  model  .  nwheels  )  ;  7«  preallocate  for  speed 
rl2  =  zeros  (12 , 12  ,  model  .  nwheels  )  ;  7#  preallocate  for  speed 
R  =  zeros  (model  .  nwheels  *12)  ;  7»  preallocate  for  speed 
for  i  =  1 : model . nwheels 

r3  ( :  ,  :  ,  i )  =  Rmaker  1  (model.  zps(:,i));  7»  r3  is  ind.  rot.  matrix 
7»blkdiagn  creates  super-duper  4x  rot.  matrix  (x  .  .  a  .  .  dx  .  .  da  .  . ) 
rl2(:,:,i)  =  blkdiagn(r3(:,:,i),4);  7»  stores  all  rl2^s  in  array 
R  (  12*  i  -  1 1  :  1 2*  i  ,  12*  i  - 1 1  :  1 2*  i  )  =  rl2(:,:,i);  7o  giant  rotation  matrix 

end 


7«7«  mi’s-  mass  matrix  inverses 

7«  compute  M  inverses  once  and  store  for  speed 
mbusi6  =  model  .  msupport  ■“ -1 ; 

mi6  =  model  .  mwheel  ■“ -1 ; 

mi 12  =  blkdiag (eye(6),  mi6); 

7«  mi  for  dummy  mass 
if  model. dummy 

mdummyi6  =  model . mdummy " -1 ; 

end 


7«7«  A  Matrix 

7«  Create  a-matrices  and  assemble  into  A 

ss.A_rb  =  blkdiagn ( [zeros (6)  eye  (6) ;  zeros(6,12)]  , model . nwheels ) ; 

7«  7«  these  are  useful  sometimes  when  the  model  breaks 
7*  size(rzl2(:,:,l)) 

7*  s  ize  (  blkdi  ag  (  Mbus  i6  ,  Mi6  )  ) 

7«  s  i  z  e  ( k ) 

7«  K  and  C  terms 

for  i  =  2 : model . nwheels 

1  K 

K  =  rl2(:,:,i)  *  blkdiag (mbusi6 ,mi6)  *  model. k  *  rl2(:,:,i)’; 
ss.A_rb(7:12,l:6)  =  ss . A_rb  (7 : 12 . 1 : 6)  +  K(l:6,l:6): 

ss.A_rb(7:12.12*i-ll:12*i-6)  =  ss . A_rb  (7 : 1 2 , 12* i - 1 1 : 12 * i -6)  +  K(l:6,7:12); 
ss . A_rb (12*i-5 : 12*i , 1 : 6)  =  ss . A_rb  (  12* i -5 : 1 2* i , 1 : 6 )  +  K(7:12,l:6); 
ss  .  A_rb  (  12*  i -5  :  12*  i  ,  12*  i  -  1 1  :  12  *  i -6)  =  ss  .  A_rb  (  1 2*  i -5  :  1 2*  i  ,  1 2*  i  -  1 1  :  1 2*  i -6  )  + 
K(7: 12 ,7: 12) ; 


7.  C 
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C  =  rl2(:,:,i)  *  blkdiag  (mbusiS  ,ini6)  ♦  model. c  *  rl2(:,:,i)’; 
ss . A_rb (7 : 12 ,7 : 12)  =  ss . A_rb  (7 : 12 , 7 : 12)  +  C(l:6,l:6); 

ss.A_rb(7:12.12*i-5:12*i)  =  ss . A_rb (7 : 12 , 1 2* i -5 : 1 2* i )  +  C(l:6,7:12); 

ss . A_rb (12*i-5 : 12*i ,7 : 12)  =  ss . A_rb  (  12* i -5 : 1 2* i , 7 : 1 2 )  +  C(7:12,l:6); 
ss  .  A_rb  (  12*  i -5  :  12*  i  ,  12*  i -5  :  12*  i  )  =  ss  .  A_rb  (  12*  i -5  :  1 2*  i  ,  12*  i -5  :  12*  i  )  +  ... 

C (7 ; 12 ,7 : 12)  ; 

end 

ss . A_rb  (7 : 12 , 1 : 6)  =  ss . A_rb  (7 : 12 , 1 : 6)  +  mbusi6  *  model. busk; 

here  is  my  dummy  mass  for  vibe 
if  model. dummy 

x=model . nwheels  *12  +  1 ; 

a  =  [zeros  (6)  eye(6)  ;  mdummyi6 *model . kdummy ( 7 : 1 2  ,  7 : 12 )  zeros(6)]; 
ss . A_rb  =  blkdiag (ss . A_rb , a) ; 

ss . A_rb (7 : 12 , 1 : 6)  =  ss . A_rb (7 : 12  ,  1  :  6)  +  mbus i6 *model . kdummy  ( 1 : 6  ,  1  :  6 ) ; 

ss . A_rb (7 : 12 , X : x+5)  =  ss . A_rb (7 : 12 , x : x+5)  +  mbus i6 *model . kdummy ( 1 : 6 , 7 : 1 2 ) ; 

ss . A_rb (x  +  6 : x  +  1 1  ,  1  :  6)  = ss . A_rb (x  +  6 : x  +  1 1  ,  1  :  6)  +  mdummy i6 *model . kdummy ( 7 : 12 , 1 : 6 )  ; 

clear  a 

input. qO  =  [input. qO;  zer os ( 1 2 , 1 ) ] ; 
end 


%%  B  Matrix 

ss.B  =  zeros  ( 12*model  .  nwheels  ,  model  .  nwheels  )  ;  “/o  preallocate  for  speed 
clear  a 

for  ii  =  2 : model . nwheels ; 
a  =  zeros (12,2); 

a([7  10], 1)  =  input .b([7  10], ii); 
a([8  11], 2)  =  input .b([8  11], ii); 

ss . B  (  1 2* ii - 1 1 : 1 2* ii , 2* ii - 1 : 2* ii )  =  rl2(:,:,ii)  *  mil2  *  a; 

end 

clear  a 

%  input. bl  creates  constant  direction  input  forces  (just  for  testing) 
input. bl  =  mil2*input . bl ; 
input . bl  =  input . bl ( : ) ; 
ss.B  =  [ss.B  input. bl]; 

7,7,  W  matrices 

7,  these  will  actually  be  created  in  the  solver  ,  but  I  do  as  much  work  here 
7,  as  possible  to  save  integration  time 

7.  W 

7,  stretch  vector:  mathworks.com/matlabcentral/newsreader/view_thread/129069 
input . w_  (  1  ,:  )  =  input . w  (  1  ,:)  ; 

input . w_  (2  ,  :  )  =  ( input . w  (2  ,  :  )  -  input . w  (  1  ,  : ) ) / ( input . time )  ; 

a  =  repmat (input.w_(l,:)  ,12,1); 
b  =  repmat (input.w_(2,:)  ,12,1); 

7,  input  .  w_  =  input  .  w_  ’  ; 

input .W(:,:,l)  =  diag(a(:)); 
input .W(:,:,2)  =  diag(b(:)); 
clear  a  b 

7,7«  0  matrices 

7,  these  will  actually  be  created  in  the  solver  ,  but  I  do  as  much  work  here 

7,  as  possible  to  save  integration  time 

Ip  =  model . mwheel (6 , 6) ; 

ss.G  =  zeros ( 12*model . nwheels ) ; 
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for  i  =  2 : model . nwheels 
j  =  12*i; 

ss.G(j-2:j,j-2:j)  =  Ip  *  r3(:,:,i)  *  mi6(4:6,4:6)  * 
[[0  -1  0];[1  0  0] ; [0  0  0]]  *  r3(:,:,i)’; 

end 


’///«  find  integration  times 

7,  Each  wheel  must  be  integrated  to  and  from  the  time  when  it  crosses 
7«  theta  =  pi/2.  Here  I  find  those  times  and  create  a  time  vector  which 
%  contains  the  appropriate  integration  limits . 

clear  a  b 

a  =  [input.w_(2,:)/2;  input.w_(l,:);  -pi/2*ones(l, model. nwheels)]^; 

b=[]  ; 

for  ii  =  1 : model . nwheels 

b  =  [b  roots(a(ii,:))^];  %  only  need  the  max  time 

end 

7«  remove  solutions  outside  of  my  time  window  (0-tf) 
b=b(b>0) ; 

b=b (b  < input .time); 

t  =  unique ( [0  b  input . t ime ]) ;  %  [0,  intermediate  times,  tf] 

clear  a  b 

7«7«  fix  matrices  for  dummy  mass 
if  model. dummy 

ss.B  =  [ss.B;  zeros ( 12 , s ize ( ss . B  ,  2 ) ) ]  ; 
ss.G  =  blkdiag ( ss . G , zeros ( 1 2) ) ; 

X  =  model . nwheel s *  1 2 ; 

input. W  =  cat (2 , input .W , zeros (x , 12 , 2) ) ; 
input. W  =  cat ( 1 , input . W , zeros ( 12 , x  +  12  ,  2 ))  ; 
clear  x 

end 

7«7«  solve  equation 

%  run  with  state_space_model  (no  need  to  run  the  entire  model  10000  times) 

options  =  model . opt  ions ; 

output  .  t  =  []  ; 

output  .  q  =  []  ; 

for  ii  =  l:length(t)  -  1 

time  =  [t(ii)  t ( ii +1 ) - 1 E-80] ; 

[  output  _t  , output _q]  =  ode45(@state_space_model,time, input . qO , opt ions  , input , ss ) ; 
output. t  =  [output. t;  output_t]; 
output. q  =  [output. q;  output_q] ; 

input .qO  =  output.q(size( output . q , 1 ) , : ) ; 
end  ; 

end 

function  dq  =  state_space_model (t , q , input , ss  ) 

7.  global  A  q  B  u 
%  global  peanut  teakettle 

%  W(t)  varies  from  wO  to  wf -  used  to  get  correct  gyro  stiffness  in  A 
W  =  input.W(:,:,l)  +  input .W(:,:,2)  *  t; 

7«  vibe  input  has  to  look  like  a  cosine  (start  from  1).  If  it  looks  like  a 
7«  sine,  the  position  drifts.  Since  the  circular  force  of  the  imbalance 
%  needs  both  sines  and  cosines  ,  I  have  to  include  sines  ,  but  I  accomplish 
%  it  by  not  starting  them  until  theta  =  pi/2 

7«  this  is  where  i  determine  theta 
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theta  =  (  input  .  w_  ( 1  ,:)  *t  +  input  .  w_  (2  ,  : ) /2*t ‘“2)  ; 

%  here  is  wheel  speed  at  t=time 
w  =  ( input . w_  (  1  ,:  )  +  input . w_  (2  ,  :  ) * t  )  ; 

7,  peanut  =  [peanut  w(2)]  ;  see  w  while  testing 

7.  teakettle  =  [teakettle  t]  ;  see  t  while  testing-  to  match  up  with  w 

7«  this  is  the  magnitude  of  the  centripetal  force  input 
u_mag  =  input  .  e  *  w.'‘2; 

%  and  here  1  test  and  only  include  the  sine  if  theta>=pi/2 

7.  sin/cos  (theta)  is  direction;  u_mag  is  magnitude 

u  =  [cos (theta) . *u_mag ;  sin ( thet a ) . * u_mag . * ( abs ( thet a ) >=pi /2) ] ; 

u  =  [u(:);  1]  ; 

B  =  ss . B  ; 

A  =  ss.A_rb  +  W*ss.G; 

dq  =  A*q  +  B*u ; 
end 


Listing  B.2:  ../models/casejreader.m 

function  varargout  =  case_reader (f ile_num) 

7*  clc  ; 

%  [model , input ]  =  case_reader ( f ile_num) 

%  show  =  case_r eader ( f ile_num ) 

% 

7«  when  called  with  2  output  arguments,  returns  model  and  input,  which  are 
%  structures  that  contain  everything  model. m  needs 

7«  when  called  with  1  output  argument,  returns  show  (wheels  and  axes) 

%  check  for  correct  number  of  arguments  (1  or  2) 

if  nargout~=l  &&  nargout  ~=2;  error (' something  is  horribly  wrong’);  end; 

filename  =  [’cases/’  num2str (f ile_num)  ’.txt’]; 

check  =  struct(’wh eels’, 0,’ ax es’,0); 
local . L  =  0 ; 


%%  reads  file  only  to  find  show-  which  wheels  and  axes  to  display 
7«  show  contains  show,  wheels  and  show,  axes 
if  nargout  ==  1 

fid  =  f open ( f i lename ) ; 
while  ''feof(fid) 

key  =  textsc  an  (fid,  ’7«  s’, 1,  ’delimiter’,’  ’); 
value  =  textsc  an  (fid,  ’7«  s’, 1,  ’delimiter’,  ’\n’); 

%  char  (key {1  , 1})  7o  for  testing 

switch  char (key{l  ,  I}) 
case  ’ wheels :  ’ 

if  check . wheels ;  disp(’wheels  multiply  defined’); 
show. wheels  =  str2num ( cell2mat (value {:})) ; 
check. wheels  =  1; 
case  ’ axes : ’ 

if  check. axes;  disp(’axes  multiply  defined’);  end 
show. axes  =  st r2num ( cell2mat ( value {:})) ; 
check . axes  =  1 ; 
otherwise 

7o  do  nothing 


end 


end 


end 

f close  (fid) ; 

varargout (1)  =  {show}; 
return 
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end 


%%  nargout  must  be  2->  read  file  for  model  configuration 
%%  set  up  geometry  of  wheels  (# ,  orientation) 

7,  orientation  of  spin  axis  (wheel  z)  in  global  coords 
model,  zps  =  [[0  0  1]’  [0  0  1]’  [0  0  -1]']; 


7«  #  wheels  derived  from  zps 
model . nwheels  =  size (model . zps  ,  2) ; 

%%  set  variables  to  defaults 

model,  busk  =  zeros  (6)  ;  7o  no  bus/support  spring 
model  .  dummy  =  0;  7,  no  dummy  mass  (appendage) 

7»  define  dummy  mass  terms  anyways,  we’ll  just  ignore  them  if  dummy=0 
model. mdummy  =  diag([l  1  1  10  10  1]); 

kdummy  =  1000*pi  “2*diag  (  [0  0  0  1  1  0]  );  7» *  1  ■  23456  ; 
model . kdummy  =  [-kdummy  kdummy;  kdummy  -kdummy]; 

input  .  w  =  zeros  (2  ,  model  .  nwheels  )  ;  7«  w  =  0  for  all  wheels 


7»  default  is  no  initial  input  or  state 
qbuilder  =  zeros ( 12 ,( model . nwheel s +model . dummy )) ; 
input . b  =  qbuilder; 
input. bl  =  qbuilder; 

7»  default  integration  time 
input . t ime  =  1 ; 


7»  empty  ode45  options 
model . opt  ions  =  [] ; 


%%  read  the  case  file  and  set  up  the  model 
fid  =  f open ( f ilename ) ; 
while  "■feof(fid) 

key  =  textsc  an  (fid,  ’7«  s’, 1,  ’delimiter’,’ 
%  char  ( key  {  :  }■ ) 

value  =  textscan  (fid  ,’ 7#s ’,  1 ,’ delimiter  ’ 
value  =  cell2mat (value {:,:}) ; 
switch  char (key {1 , 1}) 

case  ’config:’  %  complete 
7«  value{l,l}  for  testing 

7«  config  =  value; 

switch  value 

case  ’validation’ 

[k  ,  c  ,  mwheel  ,  msupport  ,  L] 
case  ’ super  ’ 

[k , c , mwheel , msupport , L] 
case  ’ thesis ’ 

[k , c , mwheel , msupport  ,  L] 
case  ’ sub  ’ 

[k , c , mwheel , msupport  ,  L] 
otherwise;  error (’ uh - oh ’ ) 


\n  ’  )  ; 

.  ’  \n  ’  )  ; 


=  valid at ion_model; 
=  super_model; 

=  super_model; 

=  sub_model ; 


end 

model. k=k;  model. c=c;  model . mwheel=mwheel ; 
model . msupport =msupport ; 

7.  local  =  set f  ield  ( local  ,’ L ’,  L )  ; 

local . L  =  L ; 
local. config  =  value; 
case  ’busk’  7o  complete 

model. busk  =  -le6*diag([0  0011  0]); 
case  ’appendage’  7o  complete 
model. dummy  =  1; 
case  ’time:’  %  complete 

input. time  =  str2num (value ) ; 
case  ’  whz  :  ’  7o  complete 

input. w  =  ones(2,3)  .♦  str2num ( value )  *  2*pi ; 

case  ’ wrpm : ’  %  complete 
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input. w  =  str2nuin  ( value )  *  2*pi  /60; 

case  ’damping: ’  %  complete 

7,  if  str2num  (  cell2mat  ( value  {:}))  =  =  0 

%  str2num ( value ) 

if  str2num ( value ) ==0 
local . c  =  0 ; 
di sp (’ damping  is  off’) 
else  dispC’i  do  not  understand’) 
end 

case  ’  input  :  ’  7#  complete 

%  this  takes  care  of  ’input. b’  and  ’input. bl’ 
eval ( [ ’ input . ’  value  ’;’]); 
case  ’  state  :  ’  7»  complete 
eval (value ) ; 

case  ’options:’  %  complete 

eval([’model. options.’  value  ’;’]); 


otherwise  7«  return  an  error  if  something  unexpected  happens 
test  =  char (key {1  , 1}) ; 
if  isempty (test ) 

elseif  strf  ind  (’ wheels  :  axes:  7«^,test) 
else 

char (key{l , 1}) 

error(’typo  in  input  file’); 

end 

end 


end 

fclose (fid) ; 

try;  cexi st = local . c ;  end 
7«  a  =  local  .  c 
%  c  =  c 

%  b  =  model . c 

7«  adjust  damping  (usually  to  turn  it  off) 

if  exist (’ cexist ’,’ var ’) ;  model . c  =  model. c  *  local. c;  end 


7«  build  qO 

input. qO  =  qbuilder(:); 

fxy  =  1;  7«  this  is  unneccessary  -  all  scaling  is  in  e  ,  so  this  is  just  1 
fab  =  fxy  *  local.!  /  4; 

local . conf ig 

if  strf ind (’ thesis  super  sub  ’, local . conf ig) 
input . b ( 10 : 1 1 ,: )  =  input . b ( 1 0 : 1 1 ,:)* f ab ; 

input  .  e  =  le-11  *  model  .  mwheel  (1  , 1)  ;  7» 
elseif  strfind(’validation’ , local. config) 
input . e  =  1 ; 

end 


varargout(l)  =  {model}; 
varargout (2)  =  {input}; 

return 


7»  imbalance  scaling  terms-  may  be  used  to  create  e  later 
7»  F_centr ipetal  =  mrw''2;  e  =  mr ;  m  =  m_wheel  ;  e  =  eccentricity 
7»  this  imbalance  causes  x/y  vibe 

7»  this  imbalance  causes  a/B  vibe 

7»  the  model  can’t  tell  the  difference,  only  scaling  is  in  input  .  b 
7»  from  page  93  of  notes 


7»  a/B  imbalance  input  should  just  be  1,  and  i  will  here  multiply  it  by 
%  fab 

end 
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7,y,  functions  to  set  up  common  wheel  properties  for  the  3  model  types 
f unct ion  [k , c , m wheel , m support ,L]  =  validation_model 

[k , msupport , L]  =  validation_model_params ; 
c  =  k/1000; 

mwheel  =  validation_mass ; 

end 

function  [k , c , mwheel , msupport , L]  =  super_model 
[k , msupport , L]  =  thesis_model_params ; 
c  =  k/500; 

mwheel  =  super_mass; 

end 

function  [k , c , mwheel , msupport , L]  =  sub_model 
[k , msupport , L]  =  thesis_model_params ; 
c  =  k/500; 
mwheel  =  sub_mass; 

end 

%%  functions  to  create  stiffness  matrices 
function  [k , msupport , L]  =  validation_model_params 
y»  distance  from  CG_support  to  CG_flywheel 
d  =  1  ; 

L  =  0.5; 

k_mag  =  2500;  */,  (stiffness  of  each  mag  bearing) 
kt  =  2*k_mag;  */,  linear  transverse  stiffness 
kz  =  0;  %  linear  axial  stiffness 

Kt  =  k_mag  *L  ■“  2/2 ;  7,  angular-  derived  from  k_mag 
Kg  =  0;  7»  about  flywheel  rotation  axis-  zero 
k_vector  =  [kt  ,  kt  ,  kz  ,  Kt  ,  Kt  ,  Kg]  '  ; 

kl  =  “diag  (k_vector ) -diag  (  [0  0  0  d~2*kt  d'‘2*kt  0]);  %  upper  left 

a  =  [0  -d*kt ;  d*kt  0] ; 

kl (1:2 ,4:5)  =  a; 

kl (4:5  ,1:2)  =  a'  ; 

k2  =  diag (k_vector )  ; 

k2  (4:  5  ,  1  :  2)  =  a; 

k3  =  k2  ’  ;  %  lower  left 

k4  =  “diag ( k_ve ct or ) ;  %  lower  right 

k  =  [kl  k2  ;  k3  k4]  ; 

7»  mass  properties  of  support  structure 
mbus  =  10; 

11  =  10; 

12  =  10; 

13  =  10; 

msupport  =  diag  ( [mbus  , mbus  , mbus  ,  II  ,  12  ,  13]  )  ; 

end 

function  [k , msupport , L]  =  thesis_model_params 
7»  distance  from  CG_support  to  CG_flywheel 
d  =  .2; 

L  =  .2; 

7o  magnetic  bearing  stiffness 

k_mag  =  1756e3;  %  N/m  (transverse  stiffness  of  each  mag  bearing) 

kz  =  0;  7»  linear  axial  stiffness-  not  modelled  here 

Kg  =  0;  7»  about  flywheel  rotation  axis-  zero 

7o  model  stiffnesses-  derived  from  mag  bearing  stiffness 

kt  =  2*k_mag;  7»  N/m  (effective  model  stiffness-  2  bear  ings  /  f  lywheel ) 

Kt  =  k_mag  *L  ■“  2/2 ;  7o  N-m/rad  (model  angular  stiffness-  from  k_mag) 

7o  stiffness  vector 

k_vector  =  [kt  ,  kt  ,  kz  ,  Kt  ,  Kt  ,  Kg]  '  ; 

7o  turn  my  bearing  stiffnesses  into  a  stiffness  matrix 
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kl  =  -diag  ( k_ve  ct  or )  -  diag  ( [0  0  0  d'‘2*kt  d''2*kt  0]);  %  upper  left 

a  =  [0  -d*kt ;  d*kt  0] ; 

kl (1:2 ,4:5)  =  a; 

kl (4:5  ,1:2)  =  a'  ; 

k2  =  diag ( k_ve ct or  )  ; 

k2  (4:  5  ,  1  :  2)  =  a; 

k3  =  k2  ’  ;  %  lower  left 

k4  =  -diag ( k_ve ct or ) ;  %  lower  right 


k  = 

[kl  k2;  k3 

k4]  ; 

y,  mass 

properties 

of  support 

structure 

mbus  = 

10;  "/,  kg 

11  =  10 

;  '/,  kg-ni‘2 

(MQI  about 

X) 

12  =  10 

;  y,  (MQI  about  Y) 

13  =  10 

;  y,  (MQI  about  Z) 

msupport  =  diag  ( [mbus  ,  mbus  ,  mbus  ,  II  ,  12 , 13 ]  )  ; 

end 

%%  functions  to  create  mass  matrices 
function  mwheel  =  validation_mass 
m  =  10; 

L  =  0.5; 
r  =  1/8; 

It  =  m/12  *  (3*r‘2  +  L‘2) ; 

Ip  =  m/2  *  r “2 ; 

mwheel  =  di ag ( [m , m , m , It , It  ,  Ip] )  ; 

end 

function  mwheel  =  super_mass 

m  =  10;  %  kg  (mass  of  flywheel) 

L  =  .2;  %  m  (length  of  flywheel) 

r  =  .15;  %  m  (radius  of  flywheel) 

It  =  m/12  +  (3*r''2  +  L“2);  %  kg-m~2  (transverse  MDI) 

Ip  =  m/2  *  r“2;  %  (polar  MQI) 

%  this  is  the  number  that  is  defined  above 
%  It  =  0.0896; 

mwheel  =  diag ( [m , m , m , It , It  ,  Ip] )  ; 

end 

function  mwheel  =  sub_mass 

m  =  10;  %  kg  (mass  of  flywheel) 

L  =  .2;  %  m  (length  of  flywheel) 

r  =  .15;  %  m  (radius  of  flywheel) 

It  =  m/12  *  (3*r'‘2  +  L“2);  %  kg-m~2  (transverse  MQI) 

Ip  =  m/2  *  r“2;  %  (polar  MQI) 

y»  for  demonstration  of  sub-sync  whirl,  redefine  It-  cases  72xx 
It  =  0.2015; 

mwheel  =  di ag ( [m , m , m , It , It , Ip] ) ; 

end 


%  model . options . MaxStep  =  2e-4; 


Listing  B.3:  .. /models/crunch. m 

function  crunch ( conf ig , varargin ) 
global  output 

if  nargin  ==  0;  config  =  ’validate^;  end; 
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if  " is numeric (config) 

fid  =  f open ( ’ case_maker . m O ; 

cases  ={} ; 

while  "feofCfid) 

garbage  =  textscan  (fid  “/oS  ^  ,  1 delimiter \n ’)  ; 
cases  =  [cases;  textscan(fid,’case  you',1)]; 

end 

f close  (fid) ; 

cases  =  cell2mat ( cases ) ; 

%  config  =  unique (b)’; 

switch  config 

case  ’ validate  ^ 

config  =  cases ( f ind ( cases < 1000 )) ; 
case  ’test  ’ 

config  =  cases ( f ind ( cases >= 1000 )) ; 
case  ’ all  ’ 

config  =  cases  ; 

end 

end 

%  return  %  for  testing 

%  if  length ( conf ig )> 1 ;  matlabpool  open  local;  end; 

%  This  code  creates  data  for  me  if  create_data  by  calling  the  model  funct 

parfor  j j = 1 : length ( conf ig ) 
ii  =  config(jj) 

[output]  =  model(ii); 

filename  =  str cat (’ data/ ’ ,  num2str ( ii ) , ’ . mat ’ ) ; 
iSave  (filename  , output ) ; 

end 

%  if  length ( conf ig )> 1 ;  matlabpool  close;  end; 
f print f ( ’ crunched ! ! ! \r ’ ) 


funct ion  iSave  ( f ilename  , output ) 
filename 

save (filename , ’ output ’ ) 

Listing  B.4:  .. /models/graph. m 

function  h  =  graph ( conf ig , show_plots , varargin) 
h=[];  global  output; 

if  nargin  ==  0;  conf ig= ’ validate ’ ;  show_plots =0 ;  end; 

%%  determine  which  cases  to  graph 
if  i snumer ic ( conf ig ) 

%  if  I  only  request  1  (or  a  few)  plots-  display  them  on  screen 
if  nargin==l||isempty(show_plots);  show_plots  =  1;  end 
else  %  read  case_maker  to  see  which  plots  to  graph 
fid  =  f open ( ’ case_maker . m ’ ) ; 
cases  ={} ; 
while  ''feof(fid) 

garbage  =  textscan  (fid  /oS 1 delimiter \n ’)  ; 
cases  =  [cases;  textscan(fid,’case  y#u’,l)]; 

end 

f close  (fid) ; 

cases  =  cell2mat ( cases ) ; 

%  config  =  unique (b)’; 

switch  config 

case  ’validate’  %  graph  only  my  validation  cases 
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config  =  cases ( f ind ( cases < 1000 )) ; 
case  ’test^  graph  only  my  test  cases 
config  =  cases ( f ind ( cases >= 1000 )) ; 
case  ’all'  %  graph  all  cases  found  in  case_maker 
config  =  cases ; 

end 

end 

save_plots  =  1 ; 

%%  crunch  if  necessary 
crunch_this  =  [] ; 

for  ii  =  config 

filename  =  [’data/’,  num2str(ii),  ’.mat’]; 
if  ~ exist ( f ilename ) ; 

crunch_this  =  [crunch_this  ii] ; 

end  ; 

end 

if  ~ isempty ( crunch_this ) ;  crunch ( crunch_this ) ;  end 

%%  create  graphs 
for  ii  =  config 

y,  load  show,  wheels  and  show,  axes 
show  =  case_reader ( ii ) ; 

filename  =  [’data/’,  num2str(ii),  ’.mat’]; 
load (filename) 

%  Create  graphs 

ncols  =  length ( show . wheels ) ; 
naxes  =  length ( show . axes ) ; 

label  =  {’$x  ~\textrm{(m)}$’,’$y  ~\textrm{(m)}$’,’$z  ~\textrm{(m)}$’,... 

’$\theta_x  ~\textrm{(rad)}$ ’  ,  ’$\theta_y  ”\textrm{(rad)}$ ’  ,  .  .  . 
’$\theta_z  '”\ t ext rm {( rad )}$’,..  . 

’$\dot{x  }  ~\textrm{(m/s)}$’,’$\dot{y  }  ”\textrm{(m/s)}$’,... 
’$\dot{z  }  ~\textrm{(m/s)}$’,... 

’$\dot{\theta_x  }  "\textrm{(rad/s)}$’,... 

’$\dot{\theta_y  }  "\textrm{(rad/s)}$’,... 
’$\dot{\theta_z}'’\textrm{(rad/s)}$’)-; 


f igure  (  i  i  )  ; 

set ( gcf  ,  ’ units  ’,  ’inches’,  ’position’, [0  0  2.4*ncols  1.4*naxes],... 
’paperpositionmode  ’  ,  .  .  . 

’ auto ’,’papersize’  ,[2.4*ncols-.5  1.4*naxes]); 
y»  this  is  pretty  obfuscated  to  me*  need  some  comments 

for  c  =  lincols  7,  this  cycles  and  creates  naxes  plots  for  each  wheel 
wheel  =  show . wheels  ( c ) ; 

for  k  =  1:  naxes  7,  this  loop  actually  creates  each  plot 
axis  =  show . axes (k) ; 

h(k,c)  =  subplot ( naxes , ncols , ncols * (k - 1 )+ c ) ; 
plot (output .t , output .q( : , axis+12* (wheel -1)) ) ; 
tf  =  output . t ( length ( output . t )) ; 
set(gca,’xlim’,[0  tf]); 

if  c  =  =  l;  ylabel  ( label  ■[  axi  s  },’ int  erpreter ’,’ latex ’)  ;  end; 
if  k==naxes ;  xlabel(’time  (s)’);  end; 

7«  put  correct  wheel  #  label  on  only  top  subplot 
if  k==l; 

%  code’s  wheel  #’s  =  reality  +  1;  (1  is  bus) 

if  wheel  ==  4 

title ( ’ Appendage  ’ )  ; 
elseif  wheel''  =  l; 
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title([’ wheel  ^  , num2s tr ( wheel)-!]) ; 
else;  title ( ^ Support ’) ; 
end 

end 


end 


end 

%  save  and/or  show  plots 

y,  create  files  of  each  plot 
if  save_plots ; 

%  MATLAB  fig  file 

saveasCgcf  ,  [fullpathC^  .\pics\figsO  ,  ’\case  ’  ,  .  .  . 

num2str  (ii)]  ,  ’figO; 
y  pdf  file  for  inclusion  in  thesis 
saveasCgcf  ,  [fullpathC'  ./pics/pdfsO  ,  ’/case  ’  ,  .  .  . 

num2str (ii)] , ’pdf’); 

y  png  file  for  quick  viewing  of  all  cases-  this  one  has  case  label 
y  suplabel ( [ ’ Case  ’ ,  num2str ( ii ) ] , ’ t ’ ) ; 

saveas (gcf , [fullpathC ’ . /pics /pngs ’),’/case’,  ... 

num2str (ii)] , ’png’); 

end  ; 

y  show_plots 

if  ~show_plots;  close;  end; 


end 

fprintfC’graphed! ! !\r’) 


Listing  B.5:  ../models/spect_me.m 

f unct ion  spect_me (con fig , f max , wheel , var argin ) 
y  this  will  make  spectrograms 
global  S  F  T  P 

if  nargin  <  3  ;  wheel  =  3;  end 

y  done  =  [7210:7212  7222  ] ; y 73 1 0 : 73 1 2] 722 1 ; 

y  if  intersect (config , done ) ; 
y  error (’ already  completed’) 

y  end 

yy  create  graphs 
for  jj  =  1 : length ( conf ig ) 
ii  =  config(jj); 

filename  =  [’data/’,  num2str(ii),  ’.mat’]; 
load (filename) 

yy 

y  wheel  =  1;  y  for  testing 

y  ii  =  1;  y  for  testing 

y  fmax  =  450;  y  for  testing 

qnum  =  1 2* ( wheel - 1 ) +4 ; 

x= output . q (:, qnum ) ;  y  alpha  of  specified  wheel  (default  2) 

samples  =  length ( output . t ) ; 
time  =  output . t ( samples ) ; 

Fs  =  samples/time; 

segments  =  30; 
poverlap  =  . 95 ; 

window  =  f loor ( s ample s / segment s ) ; 
y  nfft  =  max  (2'' nextpow2  (window  )  ,2'‘12); 
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nfft  =  2^15; 

noverlap  =  floor(window  *  poverlap)  ; 

%  figure (1)  %  for  testing 

f igure ( i i + 10000 ) ;  clf ; 

set ( gcf ,  ’ units ’inches’,  ’position’, [0  0  8  4],... 

’paperpositionmode ’  ,  .  .  . 

’ auto ’,’papersize’,[8  4]); 

%  S=spectrogr am (x, window , noverlap ,nfft ,fs) 

[S,F,T,P]=spectrogr am (x, window , noverlap ,nfft ,Fs) ; 

if  nargin  ==  1 ; 
f max  =  1000 ; 

end  ; 

Fmax  =  floor(fmax  *  1000/F  ( 1000 ) )  ; 

%  colorbar  is  messed  up  for  signals  <le-20:  so  I  multiply  by  lelO 
pcolor(T,F(l: Fmax )  ,lelO*P(l: Fmax  ,  : ) ) 

%  pcolor(T,F(l:Fmax),P(l:Fmax,:)) 

shading  flat 
colorbar 

xlabel  ’ Time  (s)  ’ 
ylabel  ’Frequency  (Hz)’ 

warning (’ Remember  to  manually  edit  the  colorbar  scale  (it"s  lelO  too  high)’) 
y//o 

%  save  and/or  show  plots 
%  create  files  of  each  plot 

%  MATLAB  fig  file 

saveasCgcf , [fullpath(’ .\pics\figs’) , ’\spec’ ,  num2str (ii)],’fig’); 

%  %  pdf  file  for  inclusion  in  thesis 

%  saveas(gcf  , [fullpath( ’ ./pics/ pdfs ’)  ,  ’/spec ’  ,  num2st r ( ii ) ]  , ’ pdf  ’ ) ; 

y,  png  file  for  quick  viewing  of  all  cases-  this  one  has  case  label 
%  suplabe 1 ( [ ’ Case  ’  ,  num2 st r ( ii ) ]  ,  ’ t ’ ) ; 

saveasCgcf , [fullpathC’ ./pics/ pngs ’ ) , ’/spec’ ,  num2str ( ii ) ] , ’ png ’ ) ; 


end 

fprintf ( ’ spec ’ ’ed!  !  !\r’) 
end 
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Appendix  C.  Inputs 

Listing  C.l:  .. /models/cases/sample. txt 

7,  sample  configuration  file 

7«  case  =  730-  this  number  is  just  a  visual  reference 

7«  all  comments  must  be  behind  7«  ^  s 

7«  do  not  leave  any  blank  rows 

config:  validation 

busk 

appendage 
stiffness :  0 
damping :  0 
time :  10 

whz:  [[0;0]  [4;12]  [4;12]] 

7«  or  wrpm : 

input:  b(10:ll,2:3)  =  le-3 
wheels :  1  2  3 
axes :  4  5 


Listing  C.2:  ../models/cases/lOS.txt 

config :  validation 
wheels :  1  2  3 
axes  :  1  2  4  5 

t ime  :  1 

whz  :  0 

state:  qbuilder (2 , 2 : 3)  =  0.01 

Listing  C.3:  ../niodels/cases/120.txt 

config:  validation 
wheels :  2  3 
axes :  410511 
t ime  :  1 

whz  :  0 

state:  qbuilder (4 , 2 : 3)  =  0.01*[1  -1] 

state:  qbuilder ( 1 1 , 2 : 3)  =  0.4*[1  -1] 


Listing  C.4:  ../niodels/cases/125.txt 

config :  validation 
wheels :  2  3 

axes  :  4  10  5  11 

busk 
t ime  :  1 

whz:  [[0;0]  [10;10]  -[10;10]] 

state:  qbuilder  (4 , 2  :  3)  =  0.01*[1  -1] 

state:  qbuilder  (  1 1  , 2 : 3)  =  0.4*[1  -1] 


Listing  C.5:  ../niodels/cases/126.txt 

config :  validation 
wheels :  2  3 

axes :  410511 
busk 
t ime  :  1 

whz:  [[0;0]  -[10;10]  [10;10]] 

state:  qbuilder (4 , 2 : 3)  =  0.01*[1  -1] 

state:  qbuilder  (  1 1  ,  2 : 3)  =  0.4*[1  -1] 


config:  validation 


Listing  C.6:  ../models/cases/lSQ.txt 
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wheels :  1  2  3 
axes  :  1  2  4  5 

whz  :  0 

t ime :  0.25 

input:  bl(10,l)  =  0.01 

options:  MaxStep=2e-4 


Listing  C.7:  .. /models/cases/700. txt 


conf ig :  the  sis 
wheels :  1  2  3 
axes :  4  5 

whz:  [[0;0]  [0;0]  [5;5]] 

input:  b(10:ll,3)  =  100 
t ime :  0.5 

options:  MaxStep=2e-4 


Listing  C.8:  ../models/cases/701. txt 


conf ig :  the  sis 
wheels :  1  2  3 
axes :  4  5 

whz:  [[0;0]  [5;5]  [0;0]] 

input:  b(10:ll,2)  =  100 

t ime :  0.5 

options:  MaxStep=2e-4 


Listing  C.9:  ../models/cases/730. txt 

7,  case  =  730 

7«  looking  for  critical  coning  speed  at  f_con  6.8  Hz 

7«  use  the  validation  model 

%  use  the  support  spring 

%  input . e  =  1 ; 

config:  validation 

busk 

wheels :  1  2  3 
axes :  4  5 
time :  10 

whz:  [[0;0]  [4;12]  [4;12]] 

input:  b(10:ll,2:3)  =  le-3 


Listing  C.IO:  ../models/cases/4001. txt 

config :  the  sis 
wheels :  1  2  3 
axes :  4  5 

wrpm:  [[0;0]  [20;20]  [20 ; 60] ] *  1 000 

time :  10 

input:  b(7:8,3)  =  1 
options:  MaxStep=2e-4 


Listing  C.ll:  ../models/cases/4002. txt 

config :  the  sis 
wheels :  1  2  3 
axes :  4  5 

wrpm:  [[0;0]  [60;60]  [20 ; 60] ] *  1 000 

time :  10 

input:  b(7:8,3)  =  1 
options:  MaxStep=2e“4 
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Listing  C.12:  .. /models/cases/4003. txt 


conf ig :  the  sis 
wheels :  1  2  3 
axes :  4  5 

wrpm:  [[0;0]  [60;20]  [20 ; 60] ] *  1 000 

time :  10 

input:  b(7:8,3)  =  1 
options:  MaxStep=2e-4 


Listing  C.13:  ../models/cases/60 12. txt 

conf ig :  the  sis 
wheels :  1  2  3 
busk 

axes :  4  5 

whz:  [[0;0]  3000/9*[l;l]  ( 3000/9  +  5 )*[  1  ;  1] ] 

t ime  :  1 

input:  b(10:ll,2:3)  =  1 
options:  MaxStep=2e-4 


Listing  C.14:  ../models/cases/60 13. txt 

conf ig :  the  sis 
wheels :  1  2  3 
axes :  4  5 

busk 

whz:  [[0;0]  3000/9*[l;l]  0*[1;1]] 

t ime  :  1 

input:  b(10:ll,2:3)  =  1 
options:  MaxStep=2e“4 


Listing  C.15:  ../models/cases/6017. txt 

conf ig :  the  sis 

appendage 

busk 

wheels :  1  2  3  4 
axes :  4  5 

whz:  [[0;0]  3000/9*[l;l]  (3000/9  +  5) *  [1 ; 1] ] 

t ime  :  1 

input:  b(10:ll,2:3)  =  1 
options:  MaxStep=2e-4 


Listing  C.16:  ../models/cases/60 18. txt 

conf ig :  the  sis 

appendage 

busk 

wheels :  1  2  3  4 
axes :  4  5 

whz:  [[0;0]  3000/9*[l;l]  400*[1;1]] 

t ime  :  1 

input:  b(10:ll,2:3)  =  1 
options:  MaxStep=2e-4 
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